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On the Simulation of Random Excitations for 
Airplane Response Investigations on Analog 
Computers 


BERNARD MAZELSKY* ann HARRY B. AMEY, JR.** 
Lockheed Aircraft Corporation 


SUMMARY 


This paper presents a method of electrical simulation which 
permits evaluations of a wide class of aircraft dynamic problems 
involving random excitations (gust loads, taxi loads, buffeting, 
fatigue, etc.), much of which cannot be treated conveniently by 
present mathematical techniques. The statistical properties 
of the random excitation—e.g., the power spectrum and prob 
ability distribution—are simulated through appropriate syn 
thesis of electrical circuitry with the output of a commercially 
available noise generator. A description of the method of syn- 
thesis is given together with the means for measuring the neces- 
sary statistical properties. In particular, the method is demon- 
strated for the simulation of non-Gaussian atmospheric turbu- 
lence and its application to the short-period response of an air 
plane over a wide range of frequencies and dampings. Non- 
linear aerodynamic restoring foree and damping characteristics 
are also considered. Finally, the following advantages over con 
ventional methods are indicated: (1) the airplane response and 
its statistical characteristics can be quickly determined for any 
variation of physical parameters, (2) the probability distribution 
of the random excitation to the airplane need not be restricted to 
Gaussian, and (3) the statistical properties of the airplane re 
sponse can be determined for nonlinear airplane characteristics 
which are not readily handled by present mathematical tech 
niques. 


SYMBOLS 


circular frequency, rad. per sec 

resonant frequency 

power spectral density 

area under power spectrum of noise gen 
erator, o;" 

power spectral density level of noise gen 
erator 

constants used in approximation of power 


spectra 


Presented at the Aeroelasticity—I Session, Twenty-Fifth 
Annual Meeting, IAS, New York, January 28-31, 1957 

*Flight Dynamics Research Group Engineer, California 
Division. 

** Dynamicist “A,” California Division 


X1, 0 


i 


amplitude and standard deviation of nois« 
generator output 

amplitude and standard deviation of fil 
tered noise 

amplitude and standard deviation of inte- 
grated filtered noise 

amplitude and standard deviation of simu 
lated atmospheric turbulence 

amplitude and standard deviation of 
simulated airplane response 

probability density and cumulative prob 
ability functions 

time, sec 

time derivative, d/dt 

damping coefficient 

rest wing force coefficient 

constant relating airplane and electrical 
frequencies 

dynamic pressure 

area of airplane wing 

mass of airplane 

mean aerodynamic chord 

radius of gyration of airplane in pitch about 
center of gravity 

density of air 

airplane forward velocity, ft./sec 

amplitude and standard deviation of gust 
velocity 

steady-state lift due to angle of attack 

steady-state pitching moment due to angle 
of attack 

steady-state pitching moment due to rate 
of change of angle of attack, 
OCy/|O(da/dt) (4/V) 

steady-state pitching moment due 
pitching velocity, OCy/[{0(d@/dt) (é/ 1 

frequency of atmospheric turbulence, 
rad./ft. 

scale of turbulence, ft 

transfer function relating gust velocity 


and airplane vertical velocity 
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Tt. = transfer function relating simulated tur- 
bulence input and simulated airplane 
response 
Zy = impedance of narrow band pass filter 
€ = nonlinear damping coefficient 
6 = nonlinear restoring force coefficient 


Subscripts 


é refers to electrical system 

a refers to airplane system 

1 refers to input 

0 refers to output or, when used with w, the 


resonant frequency 


INTRODUCTION 


_— RESPONSE of an airplane to random excitations 
such as atmospheric turbulence, buffeting, landing 
impacts, taxting over rough surfaces, etc., has been de- 
termined by assuming a stochastic process and utilizing 
the power spectral methods of generalized harmonic 
analysis for describing the statistical character of the 
Such 
analysis has been applied to several airplane response 
In all of these 


phenomenon in a concise mathematical form. 
investigations —e.g., see references 1—5. 
studies two assumptions are required: 
tem which is being subjected to a random excitation 
must be linear; and second, the description of the prob- 
ability distribution requires only the knowledge of the 
second moment or variance which can be obtained from 
If the probability 


first, the sys- 


the power spectral density function. 
distribution of the random excitation is Gaussian, then, 
for a linear system, this approach is sufficient for 
describing not only the probability distribution of the 
random excitation but also of the output response. 

For the case of linear systems subjected to non- 
Gaussian random excitations, the power spectral ap- 
proach can be used but must be extended to provide 
information on the higher order moments of the prob- 
ability distribution of the output response. This ex- 
tended approach is described in reference 6, in which it 
is shown that calculation of the higher order moments 
of the output responses requires knowledge of multi- 
variate power spectral density functions. As the order 
of the moments increases the order of variables required 
to describe the power spectral density function also 
increases. The relations between the input and output 
functions are extremely diflicult to use because of the 
complex computations required; furthermore, they are 
applicable to linear systems only. 

The limitations to linear systems and normally dis- 
tributed excitations can be easily circumvented through 
the medium of an analog computer (or, equivalently, 
an electronic differential analyzer) and a commercially 
available electronic noise generator. These electronic 
tools not only allow for the simulation of a random ex- 
citation and dynamic system with arbitrary character- 
istics but also permit the measurement of the statistical 
properties of the response in a small fraction of the time 
usually required by digital evaluation. 

The purpose of this paper is to describe a method of 
synthesis which utilizes the output of an electronic 
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noise generator to simulate a random excitation havine 
an arbitrary power spectrum and probability distriby 
tion. 
system which is either linear or nonlinear. 


This excitation can be applied to a dynami 
The elec. 
trical circuits which can be used to measure the power 
spectrum and the probability distribution are als 
described. The method is used here to simulate non 
Gaussian atmospheric turbulence (exponentially djs. 
tributed) which is in turn applied through an electron: 
differential analyzer to a simulated airplane in its short 
period response. Studies have been made over a wide 
range of airplane damping, frequency, and nonlinear 
aerodynamic characteristics. 


Basic CONCEPTS 


In this section several concepts involving the elec 
trical simulation of the statistical properties of a ran- 
dom function will be discussed. The present treatment 
is restricted to phenomena which can be described by a 
single random variable. For example, the determi- 
nation of the airplane response through atmospheric 
turbulence would be limited to defining the gust profile 
in a vertical direction only—i.e., spanwise variations o! 
the gust are the 
principles of electrical simulation are established and 


neglected. However, once basic 
understood for the one-dimensional case, the extension 
to the (or higher) case should be 
straightforward. 

The two concepts which have played the major role 


two-dimensional 


in the description of a random variable are the power 
spectrum and the probability function. The former 
concept is analogous to that of Fourier analysis of a 
nonperiodic function, say x(t), and provides some idea 
of the relative magnitudes of its various simple har- 
monic components. In order to obtain a meaningful 
description of its spectrum, when random, the Fourier 
coefficients must be treated as averaged values over a 
sufficient period of time. The modulus squared of these 
Fourier coefficients divided by the time interval, / 
is the quantity defined as the power spectrum. Phys 
ically, the power spectral density is a measure of the 
amount of energy concentrated at different frequencies 
The probability function gives the probability that the 
random variable X(f) found between the 
prescribed values, x and x + dx, denoted symbolically 


may be 


as 
Pix < x(t)<x+dx} 


The power spectrum and probability density func 
tion are related only through their respective levels 
i.e., the area under the power spectrum is the variance 
of the distribution. If the shape of the probability 
density function is known, then knowledge of the power 
spectrum uniquely determines the probability function. 
If the shape of the probability function is not known,” 
then the power spectrum does not uniquely define the 


* This usually occurs when the input to a physical system 1s 
non-Gaussian or the system is nonlinear; in both cases, the shap« 
of the output probability distribution is unknown 
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SIMULATION OF 


probability density function, indicating that additional 
properties of the random variable are required. 

The additional properties required to determine the 
probability density function can be established by ex- 
tending the concepts on which the power spectrum 
itself was established. The results of this extension, 
which are described in reference 6, involve multivariate 
power spectra from which as many moments of the 
probability density function can be determined as are 








desired. Once the desired number of moments are 
obtained, it can be assumed that the probability func- 
The 


results of reference 6 also provide a method of relating 


tion is completely defined for practical purposes. 


the moments of the probability density function of the 
random input to those of the response of a linear sys- 
tem: this likewise involves knowledge of multivariate 
power spectra. Unfortunately, the usefulness of such 
results is limited in that the relations are complex and 
dificult to evaluate; multivariate 
power spectra of random variables of interest to aircraft 


furthermore, the 
problems have not been evaluated. In contrast, the 
probability density function which can be determined 
from the multivariate power spectrum is more readily 
measurable, and some information, although limited, is 
available for random phenomena occurring in aircraft 
response problems. 

For the present, the electrical simulation of random 
phenomena appears most feasible in terms of matching 
the two functions —power spectral density and prob 
ability density. It is assumed in this paper that these 
two functions adequately describe the pertinent prop- 
erties of the random variable. The particular approach 
used herein is, of course, not unique but rather is de- 
pendent on the particular electrical equipment at the 
authors’ disposal. However, the method of approach 
is general and does establish an engineering tool which 
eliminates many of the mathematical limitations and 
permits the further extension to random variables of 
higher dimensions. 

The simulation is accomplished through the use of 
rapidly convergent successive approximations in the 
form of electrical circuits which serve to modify the 
output of a noise generator such that a random vari 
able (time history) is produced with the desired power 
spectrum and probability density. If the probability 
density of the simulated random process is normal, 
then the simulation is limited to only the power spec- 
trum since the probability density of the noise generator 
isnormal and the linear operation required to simulate 
a desired power spectrum results in a normal prob- 
ability density (central limit theorem). For this case, 
no successive approximations are required. If the 
probability density of the simulated random process 
isnot normal, then the simulation requires, in addition 
to matching the power spectrum, a nonlinear operation 
on the output of the noise generator to modify its 
normal probability density function. This nonlinear 
operation must presumably have an effect on the power 
spectrum; consequently, measurements of the power 
spectrum must be made and, if appreciably different 
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— NOISE GENERATOR SPECIFICATION 


° MEASURED VALUE 


D. OF NOILE GENERATOR 
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0 200 300 4 


ELECTRICAL FREQUENCY, & + .p.s 


FIG.| POWER SPECTRUM OF NOISE GENERATOR 


from the previous one, a correction must be made. 
This correction is in the fori of a linear operation, which 
may effect, in turn, a further alteration of the prob 
ability density function requiring an additional non 
linear operation. Thus a method of successive ap 
proximations is required; the convergence of which, 
unfortunately, depends upon the particular problem 
and choices of the approximations for simulating the 
power spectrum and probability density functions. 
Prior to a detailed discussion of this procedure, the 
individual electrical simulation of a function having an 
arbitrary power spectrum and probability density will 


be described. 


Synthesis of Power Spectrum 

The simulation of the power spectrum by electrical 
circuitry depends first on the properties of the noise 
study, the Model 101 
Its operation 


generator. In the present 


Electronic Noise Generator is used. 
depends upon a gas thyratron as a primary source of 
noise. Various filtering and stabilizing circuits are 
used to ensure an output power spectrum that is uni 
This generator was designed 
Note 
that this frequency range is too low for use with the 
The power 


form from zero to 35 eps. 
for use with electronic differential analyzers. 


direct analog computer (McCann Analog). 
spectrum as described by the specifications is plotted in 
Fig. 1 together with several measurements that were 
made by the technique described in Appendix A. 
The magnitude of the flat part of the spectrum is noted 
by the symbol B. For the purposes of the present 
study the power spectrum is approximately flat to 
sufficiently high frequencies that the step function may 
be employed to describe it mathematically as follows: 


[¢(w) /B] noise generator I (w) (1) 
where ¢(w) is the power spectral density and B is the 
level. 

Consider the problem of determining a linear im 
pedance function, Z(w), which gives rise to an arbitrary 
power spectrum when operating on the output of the 
noise generator. An elemental impedance function 
can be obtained from a resonant circuit as follows: 


Z(p) (a + bb) (d + eb + fp?) (2 
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FIG.2 MATHEMATICAL REPRESENTATION OF SINGLE 
AND DOUBLE PEAKED POWER SPECTRUMS 


This impedance, which can be represented electrically 
by either passive or active elements, is too restrictive 
for general synthesis. A more general form can be 
determined by forming a series of partial fractions as 
follows: 


> dy = b,p 


Z(p) = (3) 
N=1 dy + erp + frp* 
or by forming an equivalent product series, 
' dy + bop ‘yy Out bab + Cub? 
Z(p) = f I f f (+) 


d, + CoP - fp” . n 1 d, “Ie CnP + frp? 


The electrical counterparts of Eqs. (3) and (4) are realiz- 
able by both active and passive elements; however, one 
should note that the passive element representation 
cannot usually be successfully achieved without the aid 
of amplifiers. (These are used for isolation purposes 
to guarantee the properties of the specified elemental 
impedance functions.) The representation 
simulated by a parallel impedance network and the 


sum iS 


product representation by a series impedance network. 


— [Podmo (Pino Po) IA Pro(w w,)? 

”? [BPmo/ (Pmo-Po) | + Pro(w/)? + bm [1 —(w wy)? 
where $,, Omi, are the peak values of ¢, and w,, w;, are 
their corresponding frequencies; w,” ¢ ,, is the limit of 
w*¢ if the decay of the first peak is extrapolated to in- 
finity, and the quantity w,* ¢, is the limit of w*@ as the 
decay of the second peak is extrapolated to infinity. 
Note that ¢, — ¢,, is actually the parameter which 
As $1, — 


¢,, approaches zero the second peak becomes an im- 


describes the sharpness of the second peak. 


Ww 


9 


Ww; 
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In general, each representation has its advantages 
For mathematical 


all 
disadvantages. example, deter 
mination of the coefficients a,, },, c,, ete. is easy wit! 
the sum representation since the interrelation between| 
terms of the series is negligible, whereas the produc 
simulation may possibly require iterative procedure 
for determining the coefficients. In general, electric, 
requirements favor the product or series representa 
tion most amenable for simulation. The choi 
depends on the particular problem; for the examples 


considered in the present study the use of the elec 


as 


tronic differential analyzer is required (due to noisy 
generator and computer frequency compatibility) and 
therefore, the product series representation will be used 

If the output of the noise generator is applied to th 
impedance function given by Eq. (4) and if the power 
spectrum of the noise generator can be represented by 
Eq. (1), the resulting power spectrum will be 


a,” + 6,?w? 
d,? + (e,2 — 2d, f,)w? + fo7w' 


d(w) 


. a,* + (6,7 — 24,¢,)o* + 6,70 . 

n=14n? + (en? — Qdyrfn)w? + fp2w' : 
This representation has the advantage that the abso 
lute value of the total impedance equals the product oi 
the absolute value of the individual impedances, henc 
each peak of the power spectrum is intimately related 
to a particular elemental impedance representation 
A further desirable feature of the product representa 
tion is that for power spectrums involving a small num 
ber of terms, a close association between the coeffi 
cients and the fundamental parameters of the power 
spectrum exists. For instance, a one-peaked power 
spectrum can be expressed in the following way: 


$(w) = Gn X 


dmb. (Pm — %)| + O1(w w,)? 
[PmPo/ (Gm — Po) | + b1(w/a,)? + Gm[1 — (w/@,)?]? 


where ¢,, is the maximum value of @¢, w, 1s the corre 
sponding frequency, ¢, is the initial value of ¢, and wy"¢ 
is the limit of w*@ as w approaches infinity. A sketcl 
of this single peaked power spectrum is given in Fig. 2 
together with a pictorial representation of these param 
ters. 

A power spectrum with two well separated peaks can 


be expressed in similar manner as 


dni  (& — b1)(w/w;)4 + O12, [1 — (w/w)??? . 
* iy (OL — Pio) (w/wi)? + (wi? 7) bmi [1 — (w/@:)* FP 


Such a phenomena almost occurs in buffeting 
occurs at 


A sketch 


pulse. 
problems where a discrete energy ‘‘pulse”’ 
what is known as the shedding frequency. 
of the double peaked power spectrum is also given i 
Fig. 2 together with a pictorial representation of thest 
parameters. 

As indicated previously, physical significance of the 
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weficients in Eq. (4) for a higher number of peaked 
ower spectrums is increasingly difficult to obtain, and, 
omsequently, a mathematical fitting may require the 
ye of iteration procedures. Fortunately, for aircraft 
namic studies, input power spectrums with more than 
two peaks seem unlikely; therefore the results presented 
jerein Should suffice for general analyses. For all cases 
indicated herein, the electrical circuits which represent 
he impedance functions that permit simulation of the 
wwer spectrums are readily obtainable by standard 
rocedures which utilize active and passive elements. 
When the number of available active elements (ampli- 
fers and associated power supplies) is inadequate, 


sive elements may be substituted. 


synthesis of Arbitrary Probability Distribution 


Ia this section the nonlinear operation required to 

jmulate an arbitrary non-Gaussian probability dis- 
iibution will be discussed on the assumption that a 
random variable, VY, with normal distribution is avail- 
ible. Such a variable is obtainable from the Model 101 
Electronic Noise Generator. If successive approxi- 
mations are needed the determination of the nonlinear 
peration relating two non-Gaussian probability dis- 
tributions is necessary; however, the results of this 
ection would still be applicable with a slight modifica- 
tion. 

The usual notation for the probability functions is 
wed: p(x,) is the probability density and /(x,,) is the 
umulative probability distribution. These two func- 
tions are related as follows: 

eX, 
P, (Xn) = f - p, (é)dé P3X, < Xn} (S) 
Note that, in Eq. (8), X,, is the random variable and 
isa particular value that Y,, may or may not exceed. 

For convenience variables X, will be 


iormalized by their corresponding standard deviation, 


all random 


a 
Be) ee = x (9) 
el 
where re =| ep, (E)dé (10) 
thus, 
£5, (€)dé = 1; P,\X, < ¥ ; = 


Xr, 
J prléd— (11) 


Consider the variable .4(¢) which is normally dis- 
tributed, variance 1, and X,(¢) which possesses a pre- 


cribed non-normal distribution, and variance |. In 
rder to simulate the random function X z(t) one need 
ily apply a nonlinear operation, f, to the normally 


istributed output of the noise generator X 4(¢), as 
Xp(t) = f[Xa(0)] (12) 
where X y(t) is prescribed by 


- eX, 
P\X a(t) < XBt = f prlé)dé = Pp(xp) (15) 
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and .\,(¢) is given as 
Pi Xa(t) < xa} = 


(1/2) [1 + erf (X4/v2)] Pa(xs) (14) 


It can be shown that a unique solution for the nonlinear 
operator, /, can be obtained in terms of P, if f is re- 
stricted to being a monotonic increasing function. On 


this basis, one can write 


. =~ , 
P\ Xa < x4} = Pif(Xa) < f(xa)} (15) 


As an example, one may note that the probability that 
X" is greater than x” is the same as the probability that 
X is greater than x. 

Since the desired output 


Xp t) = f [(X4(0)] (16) 

Eq. (15) can be rewritten as 
(vc of > os . on 
P} Xa < xa} = P\ Xp < f(x4)} (17) 


If the relations (13), (14), and (17) are combined, then 


Palf(xa)) = (1/2) [1 + erf (v4/V2)] (18) 
and the solution for the operator, /, is 
f(xa) = Pa) (1/2) [1 + erf(xa/V2)]}! (19) 


where the inverse operation ?,~! is defined such that 


Pp IP R[f(xa) } = f(x) (20) 

For the case of successive approximations, Eq. (14), 
which denotes the form of P4(x4), can be replaced 
by an alternate but non-Gaussian function which is 


prescribed, say, by Pe (x4) 


Method of Successive Approximations 


As indicated previously, when a non-Gaussian prob- 
ability function is to be simulated, a nonlinear oper- 
ation is required which may affect the synthesis of the 
power spectrum. Probably the most convenient ap- 
proach for determining a solution which involves a non- 
linear operation is a method of successive approxima- 
tions. In general, the success of convergence depends 
upon the particular problem, the use of intelligent 
initial estimates, and the order or procedure of oper- 
ations. The success of the method for the problems 
dealing with airplane dynamic response apparently 
rests on the form of the power spectrum and probability 
functions which exist for such random phenomena as 
atmospheric turbulence, runway roughness, buffeting, 
etc. A review of the literature indicates that, in gen- 
eral, the probability functions associated with these 
phenomena are either Gaussian or near/y Gaussian, 
the uncertainty being primarily due to the lack of com- 
prehensive measurements. Qualitatively, the prob- 
ability functions may be non-Gaussian due to non- 
linear effects; for instance, the atmosphere experiences 
motions of large air masses which are separated by sharp 
fronts; this in turn causes a concentration of random 
phenomena which is peculiar in behavior in comparison 
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to the random phenomena of the air masses themselves. 
The addition of many such types of concentrations 
may cause the probability function of the phenomena 
to be non-Gaussian but, fortunately, not very far from 
Gaussian. The same type of reasoning appears ap- 
plicable for the runway roughness and buffeting prob- 
lems. Examples of this trend can be found in the 
literature; for instance, measured responses of highly 
damped airplanes in atmospheric turbulence have 
indicated that distribution of the loads experienced 
while near Gaussian are apparently exponential. As 
will be demonstrated in the next section of the paper, 
which deals with applications to aircraft response, 
the random disturbance which causes this type of dis- 
tribution is likewise exponentially distributed. 

The power spectral density functions associated with 
these phenomena appear to have similar forms, one 
peak of which usually occurs at a very low frequency. 
In general as the frequency increases, the level of the 
pewer spectrum decreases rapidly, usually by the in- 
verse of the square or some higher power of frequcncy. 
For atmospheric turbulence the peak in the power spec- 
trum usually occurs well below the short period mode 
of the airplane. 

In summarizing, it appears that the power spectra of 
interest are smooth functions possessing usually not 
more than one peak while the probability functions are 
approximately, but not quite, Gaussian. The non- 
linearity which is required to generate the non-Gaussian 
probability function from the Gaussian output of the 
noise generator is not very severe and consequently 
should not greatly affect the simulation of the power 
spectrum. Unfortunately, a general method of syn- 
thesis cannot be stated; however, a suggested method, 
which has been shown to be adequate for the atmos- 
pheric turbulence problem (as treated in the next 
section) is given as follows. 

(1) Assume the spectrum of the noise generator to 
be described adequately by Eq. (1). Note that the 
adequacy is ensured by the form of the power spectrum 
to be simulated—that is, ¢(w)/B ~ 1/w"; (w > 0). 

(2) Use the product series given by Eq. (5) to fit the 
desired power spectrum and determine the coefficients 
Qn, On, Cn, etc. These coefficients then establish the im- 
pedance function Z(p) given by Eq. (4) from which an 
electrical circuit can be determined for simulation. 
Apply this electrical circuit in series with the output of 
the noise generator. 

(3) Determine the nonlinear operator, f, by Eq. (19); 
this nonlinear operator can be simulated with conven- 
tional electronic devices by a procedure to be described 
in the next section of this paper. Now connect the 
circuit which performs this operation in series with the 
noise generator-impedance combination. 

(4) After the nonlinear operation, measure the re- 
sulting power spectrum and compare with the original 
spectrum. If different, compute the difference by 
taking the ratio of the true power spectral density to 
the measured power spectral density and make this 


correction with the use of Eq. (4). This correction is 
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then accomplished electrically and the resulting cirey; 


put in series with the previous circuits. 


(5) Measure the probability of the output followin, 


the correction to the power spectrum. If differen 
from the prescribed function, a correction must }y 
generated using Eq. (19), except that the quantity 
(1/2) [1 + erf (x4/+/2)] should be replaced by the 
measured probability function. This correction to th, 
probability curve must subsequently be simulated elec 
trically and the corresponding circuit is put in series 
with the preceding circuitry. 

(6) Repeat steps (4) and (5) until no corrections are 
necessary. 

This procedure appears to converge rapidly such that 
only steps (1), (2), and (3) are usually required. This 
desirable feature is due to the fact that (a) the power 
spectrum to be simulated is reasonably smooth and 
does not require complex electrical circuitry, and (h 
the nonlinear operator, /, required to produce the non 
linearity is not severe. 

If the correction factor for the power spectrum is rela 
tively flat, then the impedance function representing 
this correction could be considered as a filter with a 
damping value somewhere near the critical value 
Under these conditions, the correction factor should 
have a negligible effect on the probability function 
Conversely, if the correction to the power spectrum 
experiences a narrow peak of high amplitude, the im 
pedance function representing this correction would 
contain a filter with a damping value much smaller 
than the critical damping, and, under these conditions 
the correction would tend to change the prescribed 
probability function to Gaussian. It is for this reason 
that the order of steps (1) and (2) is suggested. 

Suggested electrical circuits for measuring the power 
spectrum and probability functions, respectively, ar 
described in Appendixes A and B. 

APPLICATION OF METHOD 

In this section the basic principles of electrical simu 
lation given previously will be applied to a one-dimen 
sional description of atmospheric turbulence. Her 
the variable of interest is the vertical velocity of th 
gust. After the power spectrum and probability func 
tions of atmospheric turbulence are prescribed in ex 
plicit form and simulated electrically, a simplifie 
mathematical representation for the short-period re 
sponse of an airplane penetrating a sharp-edged gus! 
is given and is likewise simulated electrically. Tht 
identification of electrical quantities with the airplan 
stability parameters is given in Table 1. In addition 
cases are included where the aerodynamic dampitt 
and restoring forces are nonlinear in a form which § 
compatible with actual aerodynamic properties at & 
tremely high and low angles of attack and easily ame! 
able to electrical simulation. Finally, as a result of th 
simplified circuitry, the time histories of two or mot 
airplane responses subjected to identical atmospher! 
turbulence can be recorded simultaneously, thus pe 
mitting interesting comparisons. 
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Correspondence of Turbulence and Airplane Parameters With Electrical Analog 


Turbulence and Airplane Parameters 


Wa 
L/l 
QL 
pe sD 
WAL) / Oarm 
to? L/t 
w(t) 
o; (QL)/(4/m )o%atmos. L 
(c/m)q = (gS/MI (Cr. (Cu. + Cug)/Ry(O) 
(k/m)a = (—qSé/k,? 1 Cu (p 2)SECL, vine 
Two" * (QL) /( Wen) 
bo (QL) /(Wamar)? (4L a2, 1) 
Cwel V L/(Wane:)! 
Wa/ Owa 
Wall) |v iP (25). 
(1 } 1 (p/2)S¢Cr, Cuq/M Cur, 1)°(6a/owg?) 
(1/}1 (Cu, + Cyr) Crk, )?]})°(€a/owe? 


Simulation of Atmospheric Turbulence 

This simulation is achieved by approximating the 
power spectrum and the probability distribution of 
atmospheric turbulence. For simplicity, an exponen- 
tially distributed probability function is assumed; this 
orm appears consistent with limited measurements 
and leads to a closed form solution for the nonlinear 
operator, f. Before discussing the details of its elec- 
trical simulation, consider first the problem of analyti- 
cally approximating the power spectrum. 

The power spectrum of atmospheric turbulence has 
been investigated over an extensive range of frequencies, 
and the following analytical expression has been sug- 


gested as a good approximation : 
T) = 


(1/3) + Q°L2)/f1 + 22/2 


¢(QL)/(3Le7 atmos 


(21) 


Where Q is the frequency of turbulence in radians per 


foot, L is the scale of turbulence, and atmos, iS the 
standard deviation of turbulence defined in terms of the 


power spectrum as 


Oo stmos _ f, $;(QL)dQ 


A plot of Eq. (21) is given in Fig. 3. 
investigation the electrical simulation of Eq. (21) was 
not exact, the 1/3 factor being omitted. In order to 
omit this factor and still satisfy the relationship for the 
standard deviation given by Eq. (22), Eq. (21) must 


(22) 


For the present 


be rewritten as 

(QL) (AL oz gms 7) = W2L2/[1 + WL2]? (23) 
The form of Eq. (23) permitted simple analytical checks 
on the circuitry involved in analog simulations as well 
as on output response of the linear airplane. The 
omission of the 1/3 factor represents a small distortion 
of the turbulence power spectrum at very low frequen- 
cies where it is perhaps least known and where the 
short-period response of the airplane is negligible. 


Electrical Analog 


w,/N 
N 
We 
4 
Y,(t)/o 
B 
X4(t)// 
od: (w,)/B 
(c/m), (1/N 
(k/m), (1/N? 

V (k/m). (c/2m NV Tx," (iw 
N2|(k/m),- (c/2m),?| |dolwe)/B 
o, NV k/m), c/2m).? 

V/o 
w(t) NY ‘(k/m) c/2m 
0 go 
€,/0 


Eq. (23) can be simulated electrically once several 


important electrical and physical parameters are re 


lated. For instance, the airplane circular frequency, 
w,, corresponds to the electrical circular frequency 
w,, as 


w, = w,/N (24) 


where we shall select V to be the ratio of the scale of 


turbulence to the forward velocity of the airplane 
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N=L/U (25) 
such that the parameter 22 can be represented elec- 
trically as 

QL = w, (26) 
and the electrical counterpart of Eq. (25) is 

¢;(w,)/B = w,2/[1 + w,?]? (27) 


where B is the level of the flat portion of the noise 
generator power spectrum. For purposes of compari- 
son, Eq. (27) is plotted in Fig. 3, while the consequent 
correspondence of electrical and physical parameters 
is summarized in Table 1. In all cases the standard 
deviation of the electrical parameters can be computed 
from the relation 
. 


o? = 4 (w,)de, (28) 


Comparison of Eqs. (23) and (27) yields, for the rela- 


tionship between B and 6° atmos. 
= Be it F (29) 


The correspondence of the gust velocity, w,(¢) with its 
electrical counterpart is an immediate consequence and 
is also given in Table 1. The calculation of the stand- 
ard deviation of the electrical variables in terms of the 
properties of the noise generator power spectrum is 
given in Table 2. 

The impedance or transfer function which gives 


rise to the power spectrum of Eq. (27) is 
Zo (Pe) = [De?/(1 + pe)?]>(1/Pe) (30) 


where /p, is the operator in terms of electrical time, p, = 
iw,. Note that the term p,”/(1 + p,)* actually describes 
a very low frequency ‘“‘double’’ filter and the term 1 //, 
is an integrator. The modulus squared of Eq. (30) is 
identically eq. (27). An electrical circuit which simu- 
lates Eq. (30) is given in Fig. 4 where X.(t)/o2 is the 
normalized variable resulting from the filter, and 
X;(¢) 03 is the normalized variable which results from 
the output of the integrator. 

The quantities indicated by the symbols go, gi, and 
g, in Fig. 4 are variable gains which provide convenient 
levels for measurement on a recorder. The triangular 
symbols represent amplifiers which are used for either 
generating impedance functions or for isolation pur- 
poses. The circuit given in Fig. 4 is by no means 
unique nor does it make efficient use of amplifiers. Its 


TABLE 2 


Calculation of Standard Deviations for Simulating 
Atmospheric Turbulence and Airplane Responses 


Electrical Variable Electrical Standard Deviations, a, 
X, (0) 01 = VA 
Xe (2) = VA — (37B/4) 
X; () o, = WV Bx/4 
X, (t) % =o; = VW Br + (by definition) 


( 


, / 
Y(t) V (2B/40,*) { [2 + (c/m),|/(c/m), | * 


~) 
=: 
ll 


* 


\pproximation of Eq. (54). 
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FIG.4 CIRCUIT FOR FILTERING AND INTEGRATING RANDOM Nojs¢ 


main purpose is to illustrate the close correspondenc 
that exists between the actual parameters and the ele 
trical simulation. 

Power spectral density measurements were made oj 
the variable, X;(/), using the method described jy 
Appendix A, and are plotted in Fig. 3. |The circular 
symbol indicates that X3(¢) is normally distributed 
The agreement between the measured and theoretical 
values is seen to be quite good. A typical time history 
comparison of the output of the noise generator, fi 
tered noise, integrated noise, and simulated atmospheri 
turbulence is plotted in Fig. 5. [Also noted is the dis 
tance corresponding to one electrical second and the 
relative voltage levels of the X,(¢), X.(¢), and X;( 
variables.| Comparison of the time histories clearly 
illustrates the effects of the filter and integration oper 
ations. The filtered noise exhibits little visual change 
although very low frequency variations in the noise 
are eliminated, as expected. The integrated noise 
after filtering indicates a suppression of the high fre 
quency components. The simulated atmospheric tur 
bulence curve plotted in Fig. 5 is the same as the pre- | 
ceding curve except that a nonlinear operation has been 
performed which changes the probability density func- 
tion from a normal to an exponential form, both oi 
v hich are symmetric with mean zero. The circuitry 
which performs this operation is determined as follows 

Consider that the variable X,(¢)/o4 can be generated 
from the variable X;(¢) 0; by a nonlinear operator, / 
such that the probability density of X4(t)/o4 will be 
exponentially distributed as follows: 


p(x) = (1/7 2)e idea (31) 


Note that the form of Eq. (31) was determined such 
that the standard deviation of 4/0, is one, hence the 
standard deviation, 04, corresponds to the standard 
deviation of the filtered and integrated noise, ¢ 
This correspondence is indicated in Table 2. 


> 


The cumulative distribution of Eq. (31) is given as 


P(xs) = (1/2)e¥*** forx < 0 


P(xy) = 1 — (1/2)e7 ¥*** for x > 0 


The solution for the nonlinear operator, f, on the varr 
able X; 03, can be evaluated from Eq. (19) by deter- 
mining the inverse operation P,~'| |, with the use oi 
Eq. (32). (Note that the subscript B corresponds 1 


this case to X4/ a4.) 
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FIG.9 CUMULATIVE DISTRIBUTION OF SIMULATED ATMOSPHERIC TURBULENCE 


Po] = 1/2) In2[] forx < 0 ) 


Pijal)= 


—(1/7/2)In2}1—[]} forx,>0 \ 
Substitution of Eq. (33) into Eq. (19) gives the equa- 
tion for the operator, /, on the variable X; as 


f(x3/o3) = (1/2) In [1 + erf (x3/+/2 os) | 
for x3; < O 
f(%3/03) = (—1/+2) In [1 — erf (x3/+/203) | oe 
for x3 > 0 
Note that Eq. (34) enjoys the desired property that 
St (x3 63) = =F — X3/ 03) (35 


A plot of the operator, f(X3/ 03) for X3 > 0 is given in 
Fig. 6 together with an approximation by straight line 
segments used for electrical simulation. The plot for 
X3 < 0 would be an antisymmetric image of that shown 
in Fig. 6. Also noted in this Figure are the numerical 
ordinates, abscissas, and slope changes of the straight 
line segments. These segments are simulated electri- 
cally by a set of diode limiters which conduct when at 
prescribed levels of the input variable Y;(¢), 03. A 
schematic of the diode limiter circuit which generates 
the exponentially distributed amplitudes, Y4(f) o4 is 
given in Fig. 7. The negative values of 6, b., and bd; 
account for the negative values of X(t) /o3, the slopes, 
Am, to Am,, represent amplifier gains which are achieved 
by resistor inputs and feedbacks. 

A check on the cumulative probability distributions 
of the normalized variables X,(¢) /o1, X(t) a2, X3(t) ‘as, 
and X,4(¢)/o, were made using the circuitry given in 
Appendix B. The distributions of X,(t)/o1, V.(t)+ 
o., and X;3(t)/o3; were all found to be quite normal. 
A typical comparison is shown in Fig. 8 for the case of 
the X3(t)/o3. The measured 
cumulative distribution of the normalized 
X,(t)/o4, is shown in Fig. 9 together with the exact 
The comparison is quite 


normalized variable 


variable, 


curve given by Eq. (82). 
good, indicating that the approximation of the non- 
linear operator by straight line segments is adequate 
for all practical purposes. 

The power spectrum of the variable X,4(¢) was meas- 
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ured by the method given in Appendix A to det 
mine to what extent the nonlinear operation Caus 
spectral density changes. The results of these me, 
urements are given in Fig. 3 (represented by a squy: 
symbol). The agreement of the measured values wj 
the prescribed curve is good and consequently no ¢ 
rection or iterative procedure seems warranted, 4) 
parently the degree of the nonlinearity in the operato 
f, is insufficient to cause any appreciable effect on ¢} 
power spectrum. These circumstances are difficult ; 
explain mathematically, consequently, attempts wer 
made to investigate higher order nonlinearities, name! 
a cubic, a quartic, and a quintic. Each was approx 
mated by straight line segments and the resulting pow 
spectra were measured. All of the nonlinear operation 
were adjusted such that the standard deviations of th 
output variable, X4(¢), would be identical. As the no, 
linearity became more severe the measured spectru 
tended to deviate more from the prescribed spectrun 
The deviations were such that, at low frequencies, th 
measured power spectrum was slightly below the pr 
scribed values, while, at high frequencies, the meas 
ured values were somewhat above. The area under al 
power spectra must, of course, remain the same. 
the high frequencies, differences in power spectral dens 
ties amounted to as much as 50 per cent for the quart 
and quintiec types of nonlinearity. Thus, for extre1 
cases, which are apparently not encountered in d 
scribing random inputs of interest in airplane respons 
investigations, a procedure of successive approxim 
tions may be necessary, for practical cases the pr 


cedure does not appear warranted. 


Simulation of a Linear Airplane Response 


The rigorous dynamic description of an_ airplai 
penetrating atmospheric turbulence is very complicat 
even for a rigid airplane. For this case, the equatio1 
of motion are described in detail in reference 7; sever 
calculations are also given for the two degree of fre 
dom problem (pitch and plunge). On the basis of refer 
ence 7 and many more unpublished NACA comput 
tions, an analytical approximation was developed | 
reference S for describing the acceleration at the cent 
of gravity of an airplane due to penetration of a shar, 
edged gust which is uniform in the spanwide directio1 
This approximation is presented in reference 8 and 
repeated here for convenience. 
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w(t.) = ; 9m (wh, /Zheveck): 
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where wW, is the normal acceleration of the airplane au 
Wamar 18 its maximum value, ¢, is the time in second 


(k/m), is the undamped natural frequency of the ail 
: ; ae 
plane, (c/2m), is a damping constant, and fapeax 15 U 


time required for w, to reach peak acceleration. 


should be noted that the simplified form of Eq. (» 
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10 CIRCUIT FOR SIMULATING AIRPLANE DYNAMICS 


a 


icludes the unsteady lift effects due to gust penetra- 
“on, (Kussner effect). 

Although Eq. (36) is the result of considerable sim 
sification, the circuitry required for its simulation is 
sificiently complex, due to the two part description 
iw(t,), that further simplification of Eq. (36) is war- 
ranted. The following expression retains practically 
il of the characteristics of Eq. (36) and requires little 


ireuitry : 


sin V (k/m)q — (€/2M)aq7*ta (37) 


[he only difference between Eq. (36) and Eq. (37) is in 
the time required to reach maximum acceleration. 
For airplanes with a good static margin (high frequency 
ishort-period oscillation) this time requirement should 
orrespond closely with the time elapsed in a quarter 
vele. When such differences do exist they represent 
1error in the transfer function of the airplane (fre- 
wency response) that would be noticeabie only at 
igh frequencies where the gust intensity is least. 


o:i(w.) ~ (1/e,? 


The transfer function between the airplane’s vertical 
dlocity, w, and the gust velocity, w,, can be derived 
rom Eq. (37) to give 

a. / (pb a 9 2 
; (Wa) mar (R/M)a (c/2m)q (38) 


ie ~“ (1W,_) : * 
(R/m)q + U(C/M) a®a Wa" 


Eq. (38) can be rewritten in terms of the gust variable, 
QL, and normalized as 

1... (aL) 
V (k/m),aN? — [(c/2m),N |? 


Ww mar 


l 
- : (39) 
N2(k/m), + 1QL(c/m),N — (QL?) 
ut, since QL. = w,, Eq. (39) becomes 
i *(1@,) 
MW. )marv (k/m),aN? — [(c/2m),N |? 
l 
(40) 


N2(k/m)_ + tw, (c/m)aN — w,? 


The right-hand side of Eq. (40) can be simulated very 
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simply by the following electrical transfer function: 


| 
‘ (41) 


Tx, * (iw,) . 
(k/m), + 1(c/m),a, — w,? 
where, by comparison or Eq. (40) and Eq. (41), the fol- 
lowing correspondence of electrical and physical param- 


eters is obtained: 


(c/m). = Ni(c/m)q, (k/m), = N*(k/m),¢ (42) 


Tene (eek) /Cadaas 
V (k/m), — (c 2m),? NT x,* (iw) (43) 


These relationships are also summarized in Table 1. 
The differential equation which generates the elec 
trical transfer function of Eq. (41) is 


(d>¥ dt,2) + (c/m),(dY/dt,) + (k'm),¥ = X4(t,) (40) 
and the corresponding electrical circuit which generates 
this differential equation is given in Fig. 10. A plot 
of the modulus squared of the electrical transfer func 
tion, 7 y,", is given in Fig. 11 for an undamped natural 
frequency of V (k mM), = &, = 207 at several values of 
the damping parameter (c/m),. 

A few remarks concerning the identification of the 
airplane parameters (W,) maz, (R/m)a, and (¢/m)q appear 
worthwhile. The calculation of the maximum accel 
eration in a sharp-edged gust as well as the damping 
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FIG.11 COMPARISON OF AIRPLANE TRANSFER FUNCTIONS 
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FIG.12 COMPARISON OF AIRPLANE RESPONSE FUNCTIONS 


and frequency of the short-period oscillation under the 
influence of unsteady effects is required if the response 
is to be calculated correctly. Some simplified ap- 
proach once again appears in order. Consider first 
the quantity (wW,),.,: for straight and moderately 
swept wings this parameter depends primarily on the 
mass parameter of the airplane (ratio of airplane mass 
to air mass displaced by the airplane) since the time 
to peak is so short that the airplane responds little in 
pitch. As a result (W,)m, can be approximated accu- 
rately from solutions of a one degree of freedom analysis 
involving vertical motion—e.g., see reference 9. The 
quantities of restoring force, (k/m),, and damping 
force, (C m)g, can be approximated quite well from 
conventionally determined steady-state stability de- 
rivatives, since the frequency of oscillation involved is 
sufficiently low in the short-period mode that unsteady 
effects may be considered negligible. On this basis, 
the restoring force, (k’m),, can be determined from the 


relation 
(km), = —(qgSé ky") X 

(Cua — [(p/2)S@/M] CraCu,y} (45) 
Similarly, the damping coefficient, (c m),, can be evalu- 
ated in terms of steady-state stability derivatives as 
follows: 
(c/m)q = (¢gS/MU) X 

| Cle — (Cie + Cu) (K,,)?]j (46) 


where the stability derivatives and airplane parameters 


are defined in the Symbols Section. The results , 
Eqs. (45) and (46) are summarized in Table 1, 
The properties of the output response of an airplan 


can be determined in several forms—e.g., the power 


spectrum of airplane vertical velocity due to randoy 
gust velocity, the probability of exceeding a given air 
plane vertical velocity due to a random gust velocity 
and finally a time history of the airplane vertical 
locity due to the random gust velocity. The first ca, 
be computed readily from the well known relationshjy 
bd = oi|1|" {7 
The second is extremely difficult to determine analyt; 
cally for non-Gaussian turbulence, and the third in 
volves an impractical amount of computation. Elec 
trically, the measurements of the second and third fory 
are extremely simple and rapid (especially the third 
Consider first the power spectrum of the airplan 
response. If use is made of Eqs. (23), (39), and (47 
the following form for the output spectrum of the 
airplane, ¢,(QL), is obtained: 


, (QL) O27? 
(tH, ) ?mar(4L0* atmos. /®) [1 + 2°22]? 
N2|(k'm), — (ce 2m),?] 


N*(k/m)q — (QL)?]? + [(c/m), NOL}? 


(48 


Similarly, the output spectrum determined electrically 
can be determined, by Eqs. (27), (41), and (47), to give 


@,(w,) B= fw2 [1+ a2) 


{ ) » 19 9. 9 
p1/ [w,? — w,7]? + (c/m),70,75 (49 


o 


where w,,” = (k m),. Note that Eq. (49) is plotted in 


Fig. 12 for a natural frequency w,, = 207 and for the 
three values of damping shown in Fig. 11. The rela 
tionship between the normalized output variables 


actual and electrical, is 


Wa/ Twa = Vio (OU 


This results from the fact that the probability functions 
of the actual and electrical variables are identical 
The relationship of o,,,, and o, can also be determined 
and 1s given as 


Cwq V Le /Wamar = TyNV (k/m), — (c/2m),? (al 


Eqs. (49), (50), and (51) are also summarized in Table | 
By combining the results of Eqs. (50) and (51), the 
correspondence of Y(t) and w,(¢) can be established. 


w,(t)W L/Wanne = VU)NV (km), — (c/2m),2 (52 


The calculation of ¢,, which is useful in checking the 


computer setup, is given in Table 2 in an approxinat 


form. The exact expression can be obtained by solving 
the following integral: 
i Ww, "do, 
oy = B | — 
J) [| + w,?]? 
I <9 


wo, + [(c/m),2 — 2(k/m),?| a2 + (k/m),? 
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SIMULATION OF 


fhe integral can be evaluated in closed form, and the 


jnal expression for ¢,* in terms of the airplane param- 


eters 1S 


»= !nrB 4(c/m),[(1 + o,,2)2 — (c/'m),2)?} X 


— 3) (w,,” + 1) + (c/m),?| (c/m), + 


| [eu 


2[(w..? + 1)? — a%,,2(c/m),? (54) 





For practical airplane configurations the following ap- 
roximations can be made: 

eel, wee a ice) (D9) 
Substitution of these approximations in Eq. (54) re 
which is sum 


> 


sults in a simplified expression for a,” 





marized in Table 2. 





o,2 = (wB 4w,,*) | [2 + (c m),]/(c/m),§ (56) 

The cumulative probability distribution of the air- 
plane response was measured on the computer using 
the method described in Appendix B. The results 
{these measurements are shown in Fig. 13 for an air- 
plane with a natural frequency (electrical) of 207 rad. 
sec. and three levels of damping. For convenience, 
the probability is plotted for normalized variables 
id compared with the normal and exponential prob 
ability distributions. The following important results 
were noted. 

|) When the airplane is highly damped (almost 
critically damped) the probability distribution is iden 
tical with the gust probability distribution, which, 
herein, has a prescribed exponential form. 

2) When the airplane is lightly damped the prob 
ability distribution tends towards the normal curve and, 
forexample, at a logarithmic decrement of g = —0.011, 
corresponds identically with the normal curve within 
the accuracy of the measurements. The use of the 
ogarithmic decrement is more advantageous in de- 


scribing the damping since it is a parameter which 1s 
independent of the frequency of the oscillation as well 
The formula 


and c/m 1s 


as the electrical or physical system. 
which relates the quantity, g, with k m 
given as 

g = (c/m)/V (k/m) — (c/2m)? (57) 
The observations (1) and (2) were also checked at 
107 and w, = 
This apparently 


two other electrical frequencies w,, = 
‘Qr for a similar range of dampings. 
verifies the conclusions reached in reference 6 concern- 
ing the probability distribution’s tendency toward nor- 
mal when the airplane is lightly damped. A rigorous 
mathematical justification of this result has not as 
vet been made, as evidenced by a criticism given by 
lick in reference 9. However, nonlinear components 
of atmospheric turbulence (elaborated upon by Tick) 
do not seriously affect the tendency toward a normal 
probability distribution when the airplane is lightly 
damped. 

Typical time histories of the airplane response to 
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simulated atmospheric turbulence are given in Fig. 
14 for three different values of damping. Due to the 
simplified representation of the airplane transfer func- 
iion, three airplane simulations were possible on the 
computer at the same time and thus simultaneous air- 
plane response readings could be made. Examination 
of Fig. 14 clearly illustrates the disastrous effects of 
low airplane damping during flight through atmospheric 
turbulence. Note that the V(t) 1s 
plotted to the same scale for the three levels of damping. 


magnitude of 


Simulation of a Nonlinear Airplane Response 


The type of nonlinearities to be treated herein are of 
aerodynamic origin and arise at large angles of attack, 
such aigles of attack can be experienced readily during 
passage of atmospheric turbulence. Since little analyti 
cal effort has been made on the problem of determin 
ing the response of an airplane with nonlinear stability 
parameters, a representation is presented which is 
based on the simplifications made for the linear air 
plane but with modifications to two stability derivatives 
which are known to experience strong nonlinearities. 
The parameter (w,),,,, 18 assumed to remain unchanged 
in the analysis although a correction may be applied 
through the wing lift curve slope. For the purpose of 
demonstrating these effects of nonlinear aerodynamic 
properties on the behavior of the rigid airplane the 
following aerodynamic functions will be simulated: 

C 


ta(@) = Cr,(a@ = 


« 


Cig @) Citala 


The form of the nonlinear variation of C,;, and 
Cia 1s consistent with its physical properties and leads 
once again to electrical simulation involving relatively 
simple circuitry. The strength of the nonlinearity 
depends upon the quantities «, and 6,; these two mag- 
nitudes have been normalized by a quantity o,,” which 
is the variance of the vertical velocity of the linear air 
plane. This normalization of the nonlinearity pro- 
vides values of «, and 6, which are of the order of unity. 
On the basis of Eqs. (58) and (59), the nonlinear 
damping and restoring force coefficients for the air 
plane can be determined, 
(q’S’/MU) X 


{cm Janonlinear 


1 Creal a = 0) [1 — (€g/ Gwe?) Wa” | 
(Cara + Catg)/(R"/0)? |} (60) 
(k ™ Janoniinear (—q'S'é k,”) x 
1 Cuala = 0) [1 — (6,/ep,2)@.7] — 
l(p 2) SEC aC My’ M}i (61) 


where (C/M)anoninear 20d (R/M)anontinear Must replace 


the linear coefficients given in Eq. (37). If the ex 
pression 
V(t) + (c/m),[1 — (e ¥? o,2)|¥ + 

(k/m), [1 — (6.¥?/o,?)|¥ = X4(t) (62) 
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FIG. 13 CUMULATIVE DISTRIBUTION OF SIMULATED AIRPLANE RESPONSES 
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is used for the electrical simulation of the nonlinear 
orm of Eq. (37), then the following correspondence 
exists between the physical and electrical nonlinearities: 


5a 6, 
of o\ esr . ; ‘ee eS (63) 
o 7a sae! [(p 2)SéC tuk Mq’ M( Mal} Oy” 
“4. — =-* (64) 
 . [(¢ Ma + C Mg’) Cralky €)* |; Cy 
Bgs. (63) and (64) are summarized in Table 1 together 
with the correspondence of o,,, and o, which is required 
to evaluate €, in terms of e, and 6, in terms of 6,. 

The circuit diagram which generates the nonlinear 
differential equation [Eq. (62)] is shown in Fig. 15. 
Note that two electronic multipliers are required to 
generate both nonlinearities. Unfortunately, only one 
was available at the time of the present investigation, 
and, consequently, only one nonlinearity could be 
treated at a time. Sufficient equipment was avail- 
able, however, to make direct comparisons of the linear 
and nonlinear airplanes simultaneously. This condi- 
tion is necessary when comparisons of airplane re- 
sponse to the same time history of atmospheric turbu- 
lence is desired. 

For the case of the nonlinear airplane, the output 
spectrum has no significance, since the curve would be 
afunction of the level of the input force. Two output 
functions which do have significance are the probability 
functions and the time histories. Typical cumulative 
distributions of airplane responses }(/¢) for the non- 
linear stability derivatives are shown in Fig. 16, where 
the values of damping and restoring force at zero angle 
of attack (or airplane vertical velocity) correspond to 
cm), 2 and (k/m), 207. For purposes of com 
parison, measured values for the linear airplane are also 
shown together with the two theoretical distributions, 
the normal and exponential. The output variable 


Y was normalized with respect to o,, the standard devi 


ation associated with the linear airplane. For the 
nonlinear spring, a value of 6, = 0.04 was used; for the 
nonlinear damping, a value of «, = 0.15 was used. 


These values are exceptionally low but are required in 
order to measure the cumulative distributions during an 
elapsed time which is sufficient to provide accurate 
measurements (higher values would frequently cause 
the computer to experience a divergence). A _ static 
divergence occurred when the nonlinear spring was 
used; a dynamic divergence was experienced when 
the nonlinear damping was used. Note that, for the 
example chosen in Fig. 16, the damping of the airplane 
is exceptionally low. Comparison of the results shown 
in Fig. 16 indicates that the effects of the nonlinearities 
are not severe, which is a result, primarily, of the low 
values of ¢, and 6,. The only noticeable difference is 
that the nonlinearly damped airplane exhibits a flatter 
probability curve at high amplitudes. 

It is unfortunate that higher values of the non- 
linearity could not be tolerated while the probability 


properties were determined. However, several ob- 


servations were made concerning the effect of the non- 
linear damping and restoring force. 

(1) Variations of Nonlinear Damping—In general, 
a larger value of ¢«, can be tolerated for smaller values 
of (c m),. When a large airplane response occurs and 
the damping becomes negative for an instant, the mag- 
nitude of the negative damping is low. However, as 
long as negative damping does not occur, the highly 
damped airplane is less likely to diverge and responds 
almost the same as a linear airplane. 

(2) Variations in Nonlinear Restoring Force--When 
the restoring force of the airplane is nonlinear, the 
highly damped airplane is less likely to diverge, as 
would be expected. Noticeable changes from the 
linear airplane occur only when the airplane damping 
is low. Static divergence is almost instantaneous 
when the restoring force is small. 

A typical time history comparison of the nonlinear 
and linear damped airplane is shown in Fig. 17. Also 
shown is the input corresponding to the two output 
responses. Note that, at the end of the time history, 
a dynamic divergence occurs at an input voltage which 
is not necessarily large (in comparison to other high 
input voltages where no sign of dynamic divergence 
occurs). In this case, however, the value of linear 
damping (c/m), used is extremely low. 

A typical time history comparison of airplanes with 
nonlinear and linear restoring forces is shown in Fig. 
1S for the same values of w,. and (c/m), as that used in 
Fig. 17. Near the end of the time history a static 
divergence is seen to occur for the nonlinear airplane 
at a time where a large input voltage is experienced. 
In general, the time histories of response for the non 
linear and linear airplanes are appreciably different 
when the linear damping, (c/m),, is low. For large 
values of (c/m), the airplane responses are almost 
identical except that small differences exist at high 
response values; this in turn makes the nonlinear air- 
plane more susceptible to dynamic divergence than the 
linear airplane. In general, however, nonlinearities 
will have an important effect only in airplanes having 


marginal stability. 
CONCLUDING REMARKS 


A method of electrical simulation has been described 
which permits evaluation of a wide class of aircraft 
dynamic problems involving random excitations. The 
method has been shown to possess the usual advantages 
of the analog computer—flexibility of parameter 
changes and rapidity of response measurements ——plus 
the advantage of handling nonlinear systems subjected 
to non-Gaussian random excitations which cannot be 
treated conveniently by present mathematical methods. 

As an illustration of the method, the statistical prop- 
erties of atmospheric turbulence (power spectrum and 
probability functions) were simulated electrically. 
In addition, a simplified mathematical representation 
of an airplane’s response—having freedom in pitch and 
plunge (short-period mode)—due to penetration of a 
sharp-edged gust, was given and represented electri- 
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ally. Several trends were noted which should provide 
asight for deducing the statistical properties of at- 
yospheric turbulence from measured airplane responses. 
for exainple, the probability distribution of an air- 
plane response to non-Gaussian atmospheric turbu 
lence tends toward the normal distribution for low 
and toward the input distribution for high 


jamping 
Thus, much information on the statistical 


jamping. 
roperties of atmospheric turbulence can be deduced 
om an airplane having very high damping; con 
ersely, little or no information can be obtained from 
Other observations 





irplanes having low damping. 
vere also noted for airplanes having nonlinear damping 
md restoring force coefficients. In particular, these 
studies have indicated that nonlinearities play an im 
yortant role on/y in airplanes whose dynamic stability 
margin is small. Thus, if the airplane is highly damped 


md nonlinear aerodynamic coefficients exist, their 





fect will be negligible in altering the statistical prop 
erties of the airplane’s response from those of atmos 
pheric turbulence. 

While the illustration presented here is limited, the 
principles involved are quite general and should pro 
vide a powerful tool to the aircraft dynamicist who 





must provide rational design criteria on load problems 
which are difficult to evaluate due to their random 


nature. 


\PPENDIX A — ELECTRICAL MEASUREMENT OF POWER 
SPECTRUMS 


Several methods exist for measuring the po ver spec 
trum of a random variable. For instance, the auto 
correlation of the random variable can be measured and 
its Fourier cosine transform determined. This method 
is used primarily in digital techniques. Its major 
disadvantage is the long time requirement for obtain 
ing practical results; the method invariably requires 
the use of expensive high-speed electronic digital 
computers for successful engineering application. An 
electronic method which uses a commercially available 
electronic multiplier requires much less time and has 
been found to be adequate for practical purposes. This 
method requires the passage of the random variable 
through a narrow band pass filter. Knowledge of the 
filter characteristics and the variance of the filter out 
put is sufficient to determine the spectral density of the 
random variable. First, determine the variance by 
squaring (using the electronic multiplier) and integrat- 
ing the filter output and measuring the average slope 
of the resulting time history. Second, make use of 
the following relation between the variance and the 
integral of the output spectral density: 

o,"(w, ) [ d,(w, w, )dw - 


{ ;(w) Z (iw, 1w,)|\*"dw (A-1) 
0 


lf the filter, Z,(iw, iw,) is made sufficiently narrow so 
that the spectrum of the random variable ¢,(w) can be 
assumed to be constant in the neighborhood of w,, the 


approximate equation for o,"(@,) becomes 





. 
0,-(w@,) = Q;'w,) | Z (10, 1@),,) “da (A 2?) 
7 


A filter which meets the preceding requirements and 
which provides excellent measurements of the po-ver 
spectral density function is given as 

Z (tw, to) = 1/[o,? + (ce m)w — w?] (A-3 
where ¢ m is the smallest permissible value such that 
no amplifier overloading is experienced during the 
electrical measurements of 4g,°(w,). Note that the 
filter described by Eq. (A-3) is identical to the one used 
in the simulation of the airplane transfer function given 
by Eq. (41). Consequently, the electrical circuit re 
quired to generate Z,;{iw, Iw,) is the same as the one 
given in Fig. 10. 

Substitution of Eq. (A-3) into Eq. (A-1), 
of the simplifications given in Eq. (55), leads to the 
following relation between the electrically measurcd 
standard deviation and the power spectrum of the 


making use 


random variable. 


o(w,) = (2) m)w,"*(c m)a7(e, (A-4 


By repeating the measurements of o°(w,) for several 
values of w, a complete power spectrum can be ob 
tained in a relatively short amount of time. 

For the measurement of a power spectrum of the 
form @)(w) ~ 1 w?, the theoretical error involved in 
using Eq. (A-4+) can be computed as 


Dp, exact — Poapproximate 
Doexact 


error = 
(c/m)/|2 + (c/m)] (A-5) 


Thus, if the damping required to prevent computer 
overloading during measurement of the power spectrum 
is cm 0.1, then the error in measuring the power 
spectrum of the form ¢,(#) ~ | w* is approximately 
5 per cent. By careful adjustment of the input level 
of the random variable, values of ¢/ m can be obtained 
which are within this limit. 

APPENDIX B-— ELECTRICAL MEASUREMENT OI! 
PROBABILITY DISTRIBUTIONS 


The simplest and quickest method of measuring the 
probability characteristics of a random variable is by 
measuring the cumulative distribution. This statistical 
property may be determined by using a decision circuit 
with a relay amplifier which determines when the ran- 
dom variable possesses a voltage below a prescribed 
value. Such a circuit is given in Fig. 19. A discus 
sion of its operation and use in measuring the cumula 
tive distribution follows. 

When the normalized random input, X (¢) /o, is greater 
than a selected value of x, a diode (amplifier No. 1 in 
Fig. 19) operates which impresses a voltage (20 volts is 
chosen for electrical purposes) across a relay which opens 
a summing circuit. Conversely, when X(t)/o is less 
than x, the diode fails to operate and the summing 
circuit remains closed. Since the input to the sum 
ming circuit is a fixed reference voltage, its cumulative 
or integrated output per unit time is a function of the 


(Continued on page 712) 








On 


the Vibration of Thin Cylindrical Shells 


Under Internal Pressure’ 
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SUMMARY 


The frequency spectra and vibration modes of thin-walled 
circular evlinders subjected to internal pressure are considered 
It is shown that for very thin cylinders the internal pressure 
has a significant effect on the natural vibration characteristics 
For these cylinders, particularly those having smaller length to 
diameter ratios, the mode associated with the lowest frequency 
The 


the mode associated 


is in general not the simplest mode. exact number of 


circumferential nodes, », which occur in 
with the lowest frequency, depends on the internal pressure p 
If this number » is large, it decreases rapidly with increasing p 
when p is small, and the ‘‘fundamental”’ frequency—the lowest 
frequency at each p—increases rapidly with increasing internal 
pressure. At higher values of internal pressure the frequency 
spectrum tends to be arranged in the regular manner, the fre 
quency increases with the increasing number of circumferential 
nodes, and the lowest frequency rises slowly with the internal 
pressure. 

the vibration 


modes, and structural damping of a series of thin-walled cylinders 


Experimental results on frequency spectra, 


subjected to internal pressure are briefly described. These 
results show agreement with the features predicted by Reissner’s 
“shallow-shell’”’ vibration theory. 

The effect of slight deviation of the cylinder from perfect 


circular symmetry is discussed. 


SYMBOLS 


a = radius of cylinder 

Q,, Q2, by, ete. = coefficients in Ko’, Ky’, Ko’ 

amp = relative amplitude 

A, B,C = maximum amplitude of component 
vibrations 

Oa, Bis By = strain components 

E = Young’s modulus 

f = w/2Qr = frequency, cycles per sec 

h = thickness of cylinder wall 

Ko, Ki, ete = coefficients in frequency equation 

Z, = length of cylinder 

m = number of axial half waves 

Mz, My, Mro = moments in shell, per unit length 

n = number of circumferential waves 
(number of nodes = 2) 

V., Ng, Neg, Or, Qy = stress resultants in shell, force per 


unit length 


Vs, No = part of N,, Ny, respectively, that is 
due to internal pressure 
Ny = N,/(Eh) = nondimensional axial tension due to 


internal pressure 


Received January 21, 1957. 

} The authors gratefully acknowledge the help of Carl Thiele 
and William Simms of the Jet Propulsion Laboratory, California 
Institute of Technology, for the design and installation of the 
electronic system and instrumentation. 

* Associate Professor of Aeronautics, California Institute of 
Technology, and Consultant, The Ramo-Wooldridge Corp. 

** Professor of Aeronautics, California Institute of Tech 
nology, and Consultant, The Ramo-Wooldridge Corp. 

*** Member, Technical Staff, The Ramo-Wooldridge Corp. 


650 


nr = N,/(Eh) = nondimensional circumferential tensio) 
due to internal pressure 

p = internal pressure, psig 

py = pressure in the tension-control bellows, 
psig 

t = time 

u,v, Ww = components of displacement of a point 
on the mid-surface of the shell 

Y= ae = coordinates in axial, circumferential, 
and radial directions 

kK = pw7a?/E = frequency factor; /« = f-22avV/p/E 

A = mra/L = axial wave length factor = (mean cir 
cumference) /(axial wave length 

v = Poisson’s ratio 

p = density of shell material 

¢g = angular coordinate 

w = circular frequency = 2rf 


INTRODUCTION 


F THE GREAT VARIETY of ways in which a thin 
O cylindrical shell may vibrate, only the flexural 
vibration of the wall of the cylinder, sometimes des- 
ignated as the breathing modes, will be considered. 
The cylinders are assumed to be ‘freely supported” 
in such a manner that the ends remain circular, and 
that no restraint on the axial or tangential displace- 
ment is imposed at the ends. 

The problem is not without elements of surprise 
It has been found that the natural frequencies are 
arranged in an order which has little relation to the 
complexity of the nodal pattern. But such phenome: 
non has been explained by Arnold and Warburton,' 
whose theoretical analysis and experimental investiga- 
the 
equations (or Love's “‘first approximation’’) in describ- 
On the 
other hand, the effect of internal pressure on the free 


tions demonstrate adequacy of Timoshenko’'s 


ing the vibration of unpressurized cylinders. 


vibration of thin cylinders does not appear to be 
generally appreciated. In reference 2, it is stated that 
“the results (of analysis of unpressurized cylinders) 
are also true for a cylinder subjected to uniform stress; 
thus internal fluid pressure in a container will not 
affect the certain (thicker 
cylinders) this is nearly true. But for cylinders with 


frequency.”’ In cases 
very thin walls, such as aircraft fuselage or missile 
bodies, the internal pressure has a grave effect on the 
frequency. For instance, one may easily construct 
realistic examples in the aeronautical field in which an 
internal pressure that induces a hoop stress of order 
0.2 per cent of the tensile strength of the cylinder 
causes an increase of 10 per cent in the lowest frequency 


over that of an unpressurized cylinder. The impor- 
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tance of internal pressure on the vibration of pressure 
vessels is quite evident. 

The large apparent effect of the internal pressure on 
the lowest natural frequency can be explained by the 
same fact as noted by Arnold and Warburton: that 
for cylinders with very thin walls the mode of the 
lowest frequency is in general not the simplest one. 
For a given axial wave length, the exact number of 
circumferential nodes that occur in the mode of the 
lowest frequency depends on the bending stiffness of 
the cylinder wall and the internal pressure. This 
number is in general quite large and is rather sensitive 
to the internal pressure if the cylinder is very thin 
and if the pressure is small. Thus, a relatively small 
change in internal pressure may cause considerable 
change in the mode shape at the lowest frequency. 
When the internal pressure is sufficiently large, a 
condition may be reached at which the mode at the 
lowest frequency does not change any more; then the 
fundamental frequency varies approximately as the 
square root of the pressure as is expected in analogy 
with a stretched string or membrane. 

There is a basic difference in the frequency spectrum 
of a cylinder as compared to those of some other more 
familiar examples, such as strings, membranes, beams, 
This is the difference in the spacing of 
successive frequencies. For the transverse vibration 
of a homogeneous string the successive frequencies 
For the trans- 


and plates. 


are spaced like successive integers 1. 
verse vibration of a uniform simply supported beam 
the frequencies are spaced like n? (m an integer). For 


a rectangular membrane the transverse vibration 


frequencies are spaced like 


ia ; 
m n° 
- 4 ~ \ + 

a? b? 


For a simply supported rectangular plate, they are 


(m, n integers) 


(a, b dimensions of membrane) 


spaced as 
fmn ™~ > += (m, n integers) 


Beams, membranes, and plates of other boundary 
conditions and shapes have less regular spacing than 
those mentioned above, but the general character is 
similar. In contrast to this, an examination of the 
frequency spectra of a cylinder will show that the 
successive frequencies rarely follow any simple rule. 
Often several successive frequencies are very closely 
grouped together. (A glance at Fig. 7 will suffice at 
this point.) 

The last point has an important bearing on engineer 
ing applications of the theory. Whereas knowledge 
about the lowest frequencies is often sufficient for 
beams and plates, for thin cylinders it is quite possible 
that a more complete knowledge about the entire 
frequency spectrum is essential in order that resonance 
can be avoided. 

The problem of vibration of unpressurized cylindrical 
shells was first attacked by Lord Rayleigh in 1894, 











Z 
FIG. |b — STRESS RESULTANTS 


and later by Love and Fliigge. See reviews in reference 
1. Extensive numerical results and experiments were 
first published by Arnold and Warburton, in whose 
1953 paper the effect of end conditions was discussed 
exhaustively both from the theoretical and experi 
mental points of view. Tables for frequencies and 
modes of free vibrations of unpressurized cylinders 
were published also by Baron and Bleich.* The case 
of a highly pressurized cylinder was discussed by 
Stern,‘ and by Serbin® who solved the problem for 
the lowest frequencies for “‘nearly-inextensional”’ vibra- 
tion modes by Rayleigh’s method. 

showed that a great simpli 
fication in the analysis of shell 
achieved if the tangential inertia forces can be neg 
It is then necessary to consider only one 


Recently, Reissner® ‘ 


vibrations can be 


lected. 
component of displacement the transverse, or radial, 


component. If the number of circumferential waves 
is sufficiently large, the resulting equations are pre 
cisely Marguerre’s equations for slightly curved plates. 
Reissner obtains a very simple expression for the 
frequency. 

In the present note the frequency equation will be 
derived according to the conventional shell theory, 
and the results will be compared with that of Reissner’s. 
The vibration with the lowest fre 
quency is then discussed in order to clarify the im 
The frequency 


mode associated 
portant effects of internal pressure. 
spectra are calculated and compared with experimental 
results. 

EQUATION 


THE FREQUENCY 


The basic assumptions of Love's first approximation 
will be used (see Timoshenko, reference 8, Chapter 11). 
Let a curvilinear coordinate system be chosen, as 
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shown in Fig. la. The x-axis is directed along the infinitesimal amplitude. Let V,, V,, J7,, J7., et 
generator of the cylinder, the y-axis, y = ag, is meas- be the stress resultants and couples in the cylinder 
ured clockwise in the circumferential direction, and induced by the internal pressure and/or a system oj 
the z-axis is directed inward along the normal to the static loading. For simplicity, we shall assume that 
middle surface of the shell. Let u, v, w be the com- - = 

5 : : ; iN, = const.. N., = const. 
ponents of the displacement cf a point on the mid- . 
surface of the shell in the x, y, z directions, respectively. whereas all other stress resultants and couples vanish 
The stress resultants and couples are shown in Fig. Ib. The effect of any possible bending caused by the end 
The equations of motion are obtained by adding inertia conditions will be neglected. Let u, v, w, N,, \ 
force terms — ph(0°u/Ot?), — ph(O0?v/Ot?), — ph(O°w/dt*), M,, ete., be the vibratory displacements and vibrator 
multiplied by the radius a, to the left-hand side of stress resultants and couples —i.e., the deviations from 
Eqs. (252), p. 438, of reference S. The resulting the uniform stress status induced by the static loading 
equations are nonlinear. For the vibration problem, and the internal pressure. u, v, w, N,, N,, ete., are 
it is justifiable to linearize by considering motions of assumed to be of infinitesimal amplitude. Neglecting 


small quantities of the second or higher order, we 
obtain the equations of motion: 


(ON,/Ox) + (ON,,/adg) — (N/a) [(0°v Ox0yg) — (Ow/dOx)| — ph(O7u/Or?) = 0 | 
(ON,,/Ox) + (ON,/ady) — (Q,/a) + N,(0°v/Ox*) — ph(0*v/dt?) = 0 
(00,/0x) + (00,/adg)+(N,/a) + N,(O°%w/dx?) + (.V,/a)[(Ov/adg) + (02w/adg?) | — ph(d°w/dt?) = 0| 
(OM,/Ox) + (O.\,,/ad¢g) — Q, = 0) 
—(0M,,/ox) + (OM,/ad¢g) — Q, = of nf 
Eliminating the shearing forces Q,, Q, and substituting the stress-strain relations, we obtain the basic equations 
(0°u/Ox?) + [(1 — v)/2a*|(07u/Og"?) + [(1 + v)/2a](07v/Oxdy) — (v/a)(Ow/dx) — [V1 — v*)/Eha| X 
[(0°v/Oxdg) — (Ow/dOx)] — [(1 — v?)/Eh]|ph(02u/dt?) = 0 
(1 + v)/2a]|(0?u/Oxdyg) + [(1 — v)/2](O°v/Ox*?) + (1/a?)(0°v/Oy?) — (1/a?)(Ow/Oy) + (h?2/12a?) X 
[(0%w/Ox*O¢g) + (0%w/a*?d—p*)] + (h?/12a*)[(1 — v)(0*v/Ox?) + (0°v/a*dy*)| + [(1 — v*)/Eh] X (4) 
[.V,(0°% Ox?) — ph(0v/ot?] = 0 
v(Ou/Ox) + (Ov/ad¢g) — (w/a) — (ah?/12)V4w — (h?/12)} [(2 — v)/a](d%v/dx*dy) + (dv a*dg*)} + 
fa(l — v*)/Eh]} N(0°w, Ox?) + (Ne a)[(Ov/adyg) + (0°w/adg?)| — ph(d*w/dt?)} = 0 
where V‘ denotes the operator 


[(0°/Ox*) + (0°/a70¢?) |? 


This set of equations admits solutions in the form 


u = A cos (mrx/L) cos ng cos wt . 
‘ ; e (mean circumference) 
v = Bsin (mrx/L) sin ng cos wt (5) A = mra/L = ‘al I | 
ut axial wave length 
w = C sin (mrx/L) cos ng cos wt \ —-" 


: : The coefficients K can be written as 
where .1, B, C are constants, and m, n are integers. ; 


On substituting Eqs. (5) into Eqs. (4), a set of homo- Ky’ = Ko + ait. + ah, + ayt,n. + 


geneous linear equations in A, B, C is obtained. For asn,? + ash,” 





a nontrivial solution the determinant of the coefficients K,’ = Ki + bin, + ben, + [n2r2/(1 — v?)?] X (8) 
of A, B, C in these equations must vanish. This leads nana, + [rt/(1 — ine) 
to the frequency equation Ko! = Ket [n?/(1 — v®) Ja, + [202/(1 — v2) Ja, 
2 he PO ae . 
Kw? — Ko’? Ki'x« — Ko’ = 0 6) i : i , as . 7 v 
+ Ay , ( When 7,, 7”, vanish, the coefficients Ay’, Ai’, A» 
where reduce to Ao, Ai, Ay which are given by Arnold and 
ea. ’ , Warburton* (see reference 1). The expressions are 
x = pw’a*/E = frequency factor 
/ : 9 a/ . 
V« = (freq. ine.p.s.) X 2raV p/E ae / 
*When ny = nz = 0, Eqs. (8) differ from reference 1 in a 
Let term which is caused by a difference in the coefficient of the 
Ss - O'w/Ox*d0¢ term in the second equation of (4). The effect on 
i, = N a (Eh), ft, = N,/(Eh) (4) numerical results due to this difference is negligible 
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equation can be simplified. 
Reissner® introduced the assumption that in ‘“‘pre- 
flexural’ vibrations of the cylinder the 


dominantly 
Under this 


tangential inertia force can be neglected. 
assumption the last terms in the first two equations 


of (2) are removed. The frequency equation then 


becomes a linear equation 


K, "K sad Ko’ 


0 


DV‘w + ph(0?w/dt?) = N,(0%x 


where F is a stress function which gives the membrane 
stresses due to vibration: 
Ni = —0O°F/dOxdy 


V, = OF /dy?, N, = °F/dx’, 


On linearizing (10) by omitting terms enclosed in the 
square brackets, and substituting for the transverse 
the same expression as given in Eq. (5), 
namely, 


deflection w 


and a stress function compatible with it 


w = Csin (mmx/L) cos (ny/a) cos wt \ (LD 
F = C’ sin (mrx/L) cos (ny/a) cos wt | 
a frequency equation is obtained immediately. 
k= [\4/(m? + d2)2?] + [1/121 — v?)](h/a)? X 
(n? +- X*)? + AA? + An? (12) 


Clearly, the first term on the right-hand side gives the 
influence of stretching of the shell, the second term 
the influence of the wall bending, and the last two 
terms the influence of the internal pressure. 

Eq. (12) is derived for large values of n. In the 
unpressurized case, Arnold and Warburton' give a 
simplified linear expression for the frequency param- 





THIN CYLINDRICAL 


Ox?) + V(O°w/dy"?) — (1/a)(O2F/Ox*) + [(0?w/Ox?)(O?F/Oy") + 


SHELLS 


— v)’Ko = (1/2)(1 — v)2(1 + v)At + (1/2)(1 — v)(h?/12a") X 
[((A2 + n?)4 — 2(4 v*)\4n? — Sd2n* — 2n® + 4(1 — v?)A4+ 4A2n7 + nn] 
— -°)?Ay = (1/2)(1 — v)(A? + m7)? + (1/2)(8 — v — 22*)A2? + (1/2)(1 — v)m? + (h?/12a7) X 
[((1/2)(3 — v)(A? + n?)? + 211 — v)At — (2 — v?)APm? — (1/2)(3 + v)m* + 2(1 — v)A? + 1? | 
-y)Ko = 1+ (1/2)(8 — v)(A? + n?) + (h?/12a7) [(A? + 7)? + 21 — v)dA? + n?] 
1— v)8aq, = [C1 — v)/2]n?(n? — A)? — [CA — v)/2]n4 + [pC — v)/2]M [(2 3p y*)/2|d?n? - 
vn — (h?/12a?) (OX? + n?)}n*r? + [(1 v)/2|n4 + vrn{ 4 
((1 + v)/2]rA?n?[(2 — r)A2 + n? — (A? + n?)?] (9) 
— y*)®ao = NQ(1 — v)A2 + [(1 — v)/2|[m? + (mn? — A*)?] + (h2/12a?) (A? + n?)? AZ + [CI v)/2|n?} 
— y’)8a,_ = A24[(B — v)/2)A2m2 + [CL — v)/2|n? + wry 
— y”)%ay A*AZ + [C1 — v)/2]n?74 
— y*)'a, = [(1 + v)/2]A2n7(n? — 1) 
— y)*b, = [(8 — v)/2]n4 + 202m? — n? + vd? — (h?/12a7)n7(rX? + 1?) 
— p)2bo = [(5 — v)/2]A4 + [(5 — 2v)/2|\?n? +.d? + (4?/12a7)d7(A? + 1?)? 
SIMPLIFICATION OF THE FREQUENCY EQUATION 
Under various additional assumptions the frequency 
where A,’ is somewhat different from A,’. This 


equation, however, is not identical with that given 
by Reissner’s shallow-shell vibration theory. In order 
to obtain the latter, it is necessary to specialize further 
to large values of n, so that the shell panel between 
the nodal lines can be approximated by a slightly 
curved plate. 

As mentioned before, Reissner’s shallow-shell vibra- 
tion theory*® leads to Marguerre’s equations 


(0°v/Oy") (0? F/Ox") 2(0°w/Oxdy) (0? F/OxOy) | (10) 


\/Eh)VAF = (1/a)(0°w/Ox*?) — |(0°w/Ox*)(O°w/dOy"*) — (O°w/Oxdy)*| 


(1 — v*)x that is valid for sufficiently small 


of n. 


eter A 
values 


CALCULATIONS FOR A SMALL NUMBER OF 
CIRCUMFERENTIAL WAVES 


Following Arnold and Warburton, the end conditions 
corresponding to Eqs. (5) or (11) are called ‘‘freely 
supported.” It is clear that whereas w and v vanish 
at all times at the ends, « does not, corresponding to 
the case in which no axial constraint is imposed. 
Furthermore, the radial movement w vanishes at 2n 
equally spaced lines along the circumference, but these 
lines are not at absolute rest, because the tangential 
motion does not vanish identically there. For con- 
venience, however, lines on which the radial displace- 
ment w is stationary will be called nodal lines. A 
pair of integers m, n then specifies a particular nodal 
pattern as shown in Fig. 2. 

In view of the simplicity of Reissner’s frequency 
equation, the lengthy equations given by (6) through 
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FIG 2- CIRCUMFERENTIAL WAVES 


(9) have no value except possibly for checking the 
accuracy of Eq. (12). In making such a comparison, 
note that the important parameters are the non- 
dimensional internal pressure parameters 7,, 7,; the 
frequency factor x, the circumference to axial wave- 
length ratio \, the number of circumferential waves 
n (number of nodes = 2), and the wall thickness to 
mean radius ratio /i/a. If the ends of the cylinder 
aré closed like a boiler, and no other uniform axial 


compression is imposed, then 


n, = pa/Eh, a, = pa/2Eh 


In some applications, 7, and #7, may be quite inde- 
pendent of each other. 

The accuracy of Reissner's frequency equation (12) 
may be tested first in the unpressurized case. A 
special set of calculations made for a rather thin- 
walled cylinder with a/h = 1,750 and v = 0.3, is 
given in Fig. 3, the curves labeled ‘‘exact’’ being those 
given by Eq. (6). The agreement is not bad even for 
nas small as 2. The complicated variation of the 
frequency with the axial wave length parameter i 
and the circumferential wave number 7 is typical of 
those demonstrated by Arnold) and Warburton. 
Similar curves drawn for other values of radius to 
thickness ratio a/h are presented in reference 9, in 
which it is shown that the frequency spectra are 
strongly influenced by the radius to thickness ratio 
only for small values of \, provided that a/h is large. 

Further numerical comparison has been given by 
Reissner!? for the case in which the wall bending is 
completely neglected. The accuracy of the shallow- 
shell theory is shown in reference 10 to give frequencies 
within 7 per cent of the exact values for \ between 0 
to m at m = 2, for pressurized as well as the unpres- 
surized cases. 

Direct comparison between the results given by 
Eqs. (6) through (9) and those by (12) has been made 
for several cases and presented in reference 9. Of 
course, Eq. (6) gives three roots for each nodal pattern, 
whereas Eq. (12) gives only one. Comparison can 
only be made between Reissner’s expression and the 
smallest of the three roots from the exact formula. 

The examples given in reference 9 show that the 
variation of «x with #,, #7, is small for small ”. The 
smallest root « begins to vary significantly with 7,, 
i, when n > 3, but the two larger roots do not vary 


much with ’,. The nondependence of the larger roots 
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on the internal pressure is simply explained by the jay 
that these roots correspond to essentially tangent, 
motions, for which the hoop stress has little stiffenjy 
effect. It is further noted that the two larger roots Ps 
considerably in excess of the smallest roots in a yer 
For n > 2, a/h > 100, a, < 5) 
10~%, the frequencies for the tangential modes are 4 


wide range of 71,. 
¥ 


least 60 times higher than that of the predominant) 
radial mode. . 
For modes of predominantly radial motion, the orde; 
of magnitude difference between (6) and (12) is some 
what larger but similar to those shown in Fig. 
approximation certain inert; 


Since in Reissner’s 


forces are neglected, it is expected that the derive 


frequency will be higher than the exact value. This 


is readily seen to be the case. 

In the next section we shall show that for thin-walle 
pressurized cylinders the value m corresponding to th 
lowest frequency is considerably in excess of 2. Hene 
one may conclude that Reissner’s expression is satis 
factory for engineering applications in which the lowe: 
frequencies are of major interest. 


Tue MopeE ASSOCIATED WITH THE 
LOWEST FREQUENCY 


It is of interest to observe some simple deductior 


. with regard to the lowest frequency and mode. 


According to Eq. (12) for given values of \, a) 


na, and n,, the frequency factor « reaches a minimum 
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FIG 4- VARIATION OF THE SQUARE OF THE VIBRATION FREQUENCY 
WITH INTERNAL PRESSURE ,FOR VARIOUS CIRCUMFERENTIAL 


WAVE NO n 
when 
dx /On = 2n} —[2A4/(n* + d*)5] + [a?/6(1 — v*)] X 
(n? + A2) + a,f =O (13) 
If, = 0, this gives 
1 = V 121 y*)] “AVa/h — (14) 


If, in addition, \ < /a/h, then 


V 2r(a/h) (15) 


Kmin = [A?/V 3(] v*) |(h/a) (16) 
Of course, m must be an integer and the above mini 
mum value may not necessarily be reached. 

It is interesting to note that when p = 0 the funda- 
mental mode is determined by a balance between the 
bending and stretching energy. The stretching energy 
due to membrane stresses decreases with increasing 7, 
whereas the bending energy in the wall due to change 
of curvature increases with increasing m. For very 
thin shells n has to be quite large so that the con 
tributions to the frequency due to bending and stretch- 
ing are equalized. 

When the 
variation of the 
at the lowest frequency is best illustrated by an example. 
Fig. 4 shows a special case in which the axial wave 
length parameter \ = 1, the bending stiffness param- 
eter = (1/12)(h/a)? = 10-8, vy = 0.3, and a2, = 2M,. 
It is seen that m, the circumferential wave number 


internal does not vanish the 


number of 


pressure 


circumferential waves 


corresponding to the lowest frequency, varies with the 
internal pressure parameter fg; at first very rapidly, 
then slowly, as shown in Table 1 and Fig. 5. Table 
lis based on Eq. (12); therefore the range of 7, cal- 





OF THtn C_LINDAIC£AL 


SHELLS 655 


culated for n < 3 is probably quite inaccurate. 

It is clear from (12) and (13) that the lowering of n 
for the fundamental mode is due to the increasingly 
larger contribution of ig to the energy due to change of 
curvature. 

To give some physical feeling about this example, 
it may be said that A = | corresponds to a cylinder 
whose length/radius ratio is 7 or integral multiples of 
m, and that h?/(12a?) = 10~* corresponds to a radius 
thickness ratio of order 3,000. If E = 10’, then 
nm, = 10° when the internal pressure is only 0.03 psi. 

It is clear that the larger the number of circumferen- 
tial waves, n, the faster is the rate of increase of fre 
quency with the internal pressure. The rapid increase 
in the lowest frequency when 7, is small is caused by 
the fact that n is fairly large at the lowest frequency 
if the cylinder is short and if the wall is very thin. 


EXPERIMENTS 


Experiments were performed in order to verify the 
important effect of internal pressure on the frequency 
spectrum as predicted by the theory. The details of 
these experiments are described in reference 11; only 
a brief outline will be given here. 


Test Specimens 
The models (see Fig. 6) used for the tests were 
cylinders with a 3.5-in. inside diameter, made of 
24S-H aluminum alloy sheets (condenser foil) of three 
0.001, 0.002, 0.003 in., and three 
11, 7, and 3.5 in. The cylinders were 


thicknesses and 


axial lengths 
made by cementing the sheets together with lap 


TABLE 1 


Range of Internal Pressure 


n, No. of Circumferential 
Parameter, ty X 10° 


Waves at Lowest Freq 
10 0 0.1 


q 0.1 0.38 
8 0.38 0.97 
7 0.97 2.5 
) 2.5 7.0 
5 7.0 22.0 
j 22.0 94.0 
3 94.0 610.0 
y 610.0 7,000.0 
l >7 O00 
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FIG.6- SCHEMATIC DIAGRAM OF TEST SET-UP 


joints in the axial direction. The two ends were 
cemented to rings which provided a seal for internal 
pressure. 

The end rings were made of brass and were much 
heavier than the cylinders, thus ensuring circular 
cross sections at the ends. The end rings were free 
to move longitudinally and, to a certain degree, were 
free to rotate as a rigid body. 

The front end of the cylinder was sealed by a solid 
nose. The end rings of the test specimen rested on a 
2-in. diameter central shaft which had holes conveying 
the controlled internal pressure into the cylinder. 
Nitrogen bottles were used as a pressure source. It was 
possible to control the internal pressure to within 
limits of +0.01 psi. 

At the rear end of the test specinien, the end ring 
was fastened to a floating ring which was sealed against 
the shaft by a flexible bellows. To this floating ring 
were attached four bellows which were capable of 
exerting a longitudinal compressive load on the cylin- 
der, and two bellows capable of exerting a tensile load 
on the cylinder. By controlling the pressures in these 
bellows, the axial stress in the cylindrical shell was 


controlled over a wide range. 


Test Setup 

Fig. 6 shows a schematic diagram of the test setup. 
The motion of the cylinder was excited by a sound 
generator through a loud speaker. The output from 
the speaker was so adjusted as to maintain a constant 
sound level of 70 db. at a microphone situated about 
2 ft. from the cylinder. The speaker was placed about 
2 ft. on the opposite side. (Changing the speaker 
location and direction did not seem to have any 
significant effect on the frequency response.) During 
a resonant oscillation, the cylinder itself generated a 
considerable amount of sound. Since only the re- 
sultant of the sound generated by the speaker and 
the cylinder was kept constant, the excitation near 
the resonance condition was less than 70 db. In other 
words, the resonance peaks in the frequency response 
curves of the cylinder oscillation was artificially re- 


duced in amplitude. The 70 db. sound level wa 
maintained between 100 to SOOO eps. 

The motion of the cylinder was measured by th, 
variation of the gap between the cylinder and a number 
of brass buttons which rested on an inner tube 
polystyrene. The cylinder and the buttons wer 
charged and formed individual condensers, the vari 
tion in capacitance of which was measured. Th 
gap between the cylinder and the button was nominal) 
0.040 in. under no-load condition; however, this spa 
ing could not be controlled very accurately. 

For the 11-in. models, there were 24 buttons, wit} 
10 buttons equally spaced along an axial line and ty 
circles of S buttons each uniformly spaced. For tly 
7-in. models, only 12 buttons were used; for the 3.5-in 


models, only 10 buttons. 


Instrumentation and Test Procedure 


After passing through an amp ifier, the signal fron 
a button was observed in three ways: 


(1) Oscilloscope This showed the wave form 
relative amplitude, and Lissajou curves between 
signals coming from any two buttons. — Lissajou 
curves were recorded photographically and were used 
to identify the vibration modes. 

(2) Harmonic Wave .\nalyzer--A Hewlett-Packard 


_ Model 300A Harmonic Wave Analyzer was_ used 


Frequency was read from a Berkeley Universal Counter 
and Timer, Model 5510. 

(3) Voltmeter--It was found that a direct correlation 
existed between the voltmeter readings and the har 


monic wave analyzer readings. 


Each series of tests began with the unpressurized 
case. The frequency spectrum was first determined, 
the phase relationship and mode shapes were then 
studied. Then an internal pressure was selected and 
maintained constant. When the exciting frequency 
was varied, the output from the pick-up buttons 
varied. In most cases it was quite easy to pick out 
a resonance peak either on the oscilloscope or on. the 
voltmeter. But there were cases in which large re 
sponses occurred over a fairly wide band of frequencies, 
and it was somewhat difficult to decide which wer 
the natural frequencies. Occasionally, beating oc- 
curred, as could be seen on the oscilloscope. There 
were also cases in which the wave form appeared 
complicated. These difficulties were disturbing at 
first sight, but after some reflection they appeared to 
be just what should be expected, when the general 
behavior of cylinder vibration and the effect of slight 
deviation of the models from circular symmetry 1s 


considered. 
TEST RESULTS 


Frequency 
Table 2 shows a set of typical records of frequency 
spectra. In this table, ‘‘f’ denotes frequency 1 
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TABLE 2 
Frequency Spectra, Model 11-001 
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frequency, cps 
amp = relative amplitude, reading from voltmeter 
eveles per sec., ‘amp’ denotes relative amplitude in 


millivolts read from the voltmeter or harmonic analyzer, 
nd ‘‘p”’ psig—t.e., the 


pressure differential across the shell. 


is the internal pressure in 
The frequencies 
so) observed are plotted as small circles in Fig. 7. 
Fig. S shows a continuous plot of the frequency 
response of cylinder 7-001 (length = 7 in., wall thick 
ness = 0.001 in. 
with the axial-tension control bellows inserted, to a 
sund field whose resultant output was 70 db. This 


Figure is a simplification of the actual picture, which 


at a nominal pressure of 0.50 psig, 


contains many more little peaks and valleys in the 
irequeney In taking the readings, 
these small bumps were neglected. 

Each spectrum in Table 2 terminated at a reading 


response curves. 


bevond which the response signals were rather small, 
under which condition the determination of resonance 
became dubious. 

Spectra for other models can be found in reference 11. 


Mode Shape 


Altogether, modes for more than 40 points were 
identified for the 11-001 model. Fig. 9 shows an 
the approximate determination of the 
The film 


various buttons successively. 


example of 
Lissajou 
curves, from The x 
component was from button No. 3. The y-component 
was from another button (whose number is indicated 
A line in the second quadrant means in 


vibration modes. recorded the 


in Fig, 9). 
phase with button 3; that in the first quadrant means 
out of phase. The circuit for button 8 was open during 
the recording, hence the trace for that button appeared 
as a horizontal line. The attenuation for buttons 
4-24 was 10 times smaller than that used for buttons 


As mentioned previously, the nominal distance of 
(040 in. between the shell and the buttons was not 
controlled. Hence, the relative 
amplitude was not high, and plots like Fig. 9 have 


accuracy of the 


oly a qualitative significance. 
Because of the complicated mode shapes (large 1) 
lor very thin cylinders at low pressure (p = 0 or 0.05 
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psi), the determination of the modes was difficult with 
the relatively small number (24) of sensing elements 
internal however, the 


used. At higher 


modes corresponding to the lower frequencies were 


pressure, 
readily recognized. 


Damping 


From measurements of (1) the sharpness of the fre- 
quency output of the sound generator at specific 
settings, (2) the sharpness of the response of the 
cylinder at a specific forcing frequency, and (3) the 
frequency response of the cylinder—-i.e., the variation 
of peak response values with the forcing frequencies 
the damping characteristics of the model is inferred. 
The most striking feature of all these measurements 
is that the half-band width in frequency in all three 
cases did not seem to vary much with frequency range 
pressure; nor did it vary much with 
axial-tension-control 


and internal 
insertion or removal of the 
bellows. 

To interpret this 
(see reference 11) of the expected half-band width as 


result, an analysis was made 
a function of the sharpness of the forcing function, 
the damping characteristics (transfer function) of the 
cylinder, and the characteristics of the 
harmonic analyzer. It is concluded that the half- 
band width of the transfer function is about equal to 
that of the The “‘structural”’ 


damping coefficient, g, as ordinarily used in flutter 


resonance 


harmonic analyzer. 


analysis, does not remain constant; instead, the 
product, g/, does remain nearly constant for a wide 
range of frequencies. The experimental results may 


be summarized as 
gf = 6 cycles per sec. 
where / is the frequency in cycles per sec. Since g/ 


is ordinarily identified as the “‘viscous’ damping 
coeflicient, one may say that the damping characteris 
the model were of a More 


represents the amplitude of the mth 


tics of viscous nature. 
precisely, if o, 
normal mode of the cylinder, the Lagrange equation 
of motion for free oscillation may be written in the 


form 
by + Cady + @n7O, 0) (c, = 2af,g) (17) 


where w, is the natural frequency corresponding to 
¢,. Our experimental results seem to show that c, 


is a constant for all values of n. 


DISCUSSION 


The uncertainties of the experimental results came 


control. In most runs 


mainly from the pressure 
(each spectrum required about two hours), the internal 
pressure fluctuated within limits of +0.01 psig. Re- 
peated tests on different days of testing for nominally 
the same test conditions revealed that an uncertainty 
of approximately | per cent in frequency should be 


allowed. 
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Theoretical results for a ‘‘freely supported” cylinder, ends, and these frequencies can be estimated ay Tobias 
from Eq. (12), are plotted in Fig. 7. The solid curves proximately by replacing m by (m + c) in the formy,| asymmet 
refer to m = 1, the dotted curves to m = 2, and the for the freely supported ends. The constant, ¢, Jie} frequen’ 
chain dots to m = 3. Frequencies for m > 3 are between 0.2 to 0.4. The theoretical effect of the ey} cylinder, 
not plotted here. The bending modes, = 1, are conditions on the frequency in the pressurized cag dynamics 
also omitted from these plots, since in bending modes has not yet been determined. It seems evident thy] dectrom: 
the masses attached to the ends of the cylinder con- the theoretical values for the frequencies of the tes} occur un 
tribute significantly to the frequencies; whereas, for specimens should be higher than the freely supportej| around a 
n > 2, the end masses, if rigidly connected, are of small curves shown in Fig. 7, the error being largest jp] usually ‘ 
influence. The actual end conditions of the test speci- m = 1. On the other hand, the apparent mass of gj} amplitud 
mens were somewhere between freely supported and was neglected in the theoretical curves. If it wen} at reson: 
clamped. Reference 2 shows that, for an unpressur- included, the frequencies would be lowered. They} The freq 
ized cylinder, p = 0, the frequencies for the clamped factors of opposing influences have not vet bee} less diffe 
ends are higher than those for the freely supported evaluated theoretically. from dyt 
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Tobias, in reference 12, shows that the effect of small 
asymmetry in the cylinder is to resolve a single natural 
irequency into two nearby values. If a given steel 
cylinder, which in engineering practice is always 
dynamically asymmetrical, is excited by a_ small 
dectromagnetic excitor, the peak does not necessarily 
occur under the excitor. When the excitor is rotated 
around a circumference of the cylinder, two locations, 
usually 90° apart, can be determined at which the 
amplitude is the largest and at which one single peak 
at resonance is observed at each of these locations. 
The frequencies at these two locations are more or 


less different depending on the degree of deviations 
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fom dynamic symmetry. For an excitor located at 


wv other position, two peaks appear in the frequency | 


| 
Ii A re | 


the asymmetry is large, these peaks move sufficiently Hh [9 y 
, I ; he, 


== PP 


response curve, usually one larger than the other. 


ypart to appear as two distinct peaks. iif 

The explanation of Tobias’ effect is quite simple. 

We note first that for a perfectly symmetric circular 
evlinder, a set of solutions of the form 

u = .1 cos (mrx/L) sin ney COs wrt 

= Bsin (marx/L) cos ng cos wl 





frequency f (cps) 
= % 


= C sin (mmx/L) sin ng cos wt 


isadmissible, as well as that given by Eq. (5) and (11 
in fact they lead to identical frequency equations. y | | | | | | | 
Now suppose that the cylinder has a slight deviation 
irom symmetry, say, due to the addition of a small 
mass at one specific point. Then the modes (5) and 
IS) at once become associated with slightly different 


300 + + + 1 4 4+ a 


—_—— m=! 


-—— mee Theory 


irequencies. This is clear from considerations of —-—m=3 


kinetic and potential energies. Hence, in essence, the © Experiment L 


ron Y ¢ lio , ney > 0 ww » ave 
reason for a slightly asymmetric cylinder to have sode! 11-001 
resolved each frequency of a perfect cylinder into two 
is that the perfect cylinder is a degenerate case. In 200 | | | ; ; | | 
the perfect cylinder the locations of circumferential 


nodes are unspecified, for an asymmetric cylinder, 





Let the angle g in (18), which is meas- j | | | i j j 
Oo 0. 0.2 0.3 0.4 0.5 06 O7 


internal pressure p (psi) 


this is not so. 


ured from an arbitrary polar axis, be replaced by 





Substituting the mode so defined into expres 








e=~ da. 
| sons of kinetic and potential energy of a slightly FIG. 7b- MAGNIFIGATION OF THE LOWER LEFT CORNER OF FIG. 7a 
‘symmetric cylinder, we can derive an approximate NATURAL FREQUENCY VS. INTERNAL PRESSURE 
irequency which is a function of a. When a varies in 
0), 27), the corresponding frequency as a function of a 
will have a maximum and a minimum. These extrema 
correspond to the normal modes of the slightly asym | 
metric cylinder and the corresponding values of a de | | as | [| 
line the nodal line locations. Pas {ih 
The effect of slight asymmetry makes an experi- 3 | LW 
mental determination of the natural frequency by e | ' | | 
resonance somewhat difficult. Some of the extraneous : | | ength 
lrequencies that appear in Fig. 7 are undoubtedly due : ie ne 
to this origin. | | ‘| [ | | i 
The experimental results certainly showed all the 4 t { S| ‘| 
important features predicted by the theory. The tH th A nit 4 1 
discrepancy between the theory and experiment arises iy |e ii HK NT? 
irom several sources. The most important is probably aad a S00 OO oo 
the end conditions. The error caused by the shallow- s-seb tas ee ami a ak atti iii 
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shell approximation could be removed on the basis of 
existing theory, but the laborious computation has 
not been The asymmetry is 
probably responsible for some of the extraneous re- 


attempted. dynamic 
corded frequencies, and some of the deviations between 
the experimental and theoretical values. Some of the 
other experimental points that do not fall on any 
theoretical curves may perhaps be explained by higher 
values of m, because curves for m > 3 are not plotted 
in these Figures. 


CONCLUSIONS 


It appears from the foregoing that the experiments 


and the theory are in reasonable agreement. For 
cylinders with freely supported ends, Reissner’s 
approximate frequency expression is adequate for 


prediction at high internal pressures for all values of 7, 
and at lower internal pressures for nm > 3. The effect 
of other end conditions on the frequency should be a 
worth-while subject for further research. 

The experiments reveal the complicated nature of 
cylinder vibration and point to various difficulties in 
obtaining accurate data. Without a theoretical back- 
ground, it is virtually impossible to understand the 
mass of frequency readings. Furthermore, whether 
a mode will be strongly excited or not in a resonance 
test depends entirely on the particular excitation force 
field; consequently, important modes may be over- 
looked in a specific experiment. This points to the 
extraordinary importance of theoretical calculations in 
the cylinder vibration problem. On the other hand, 
the experimental results presented in this report show 
clearly that the magnitude of the response has no simple 
relation to the order of the frequency in the spectrum. 
Thus a very large response may occur at a frequency 
which is several times the lowest one and which is of 
a large order in the spectrum. Hence, in engineering 
applications of the cylinder vibration theory, it is 
imperative that the calculations be not stopped at the 
first few lowest frequencies, but that the whole spectrum 
in the frequency range of interest be examined. 

The experiments were designed for the verification 
of the main effects of the internal pressure as related 
to aeronautical engineering applications. The ac- 
curacy is not high enough to discriminate various 
refined theories of thin shells. For the latter purpose 
the reader is referred to works of Naghdi and Berry 
(1954), Kennard (1953), and Yu (1955) in recent issues 
of the Journal of Applied Mechanics. 
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SUMMARY 


‘he flow over a two-dimensional finite wedge at zero angle of 
ittack for subsonic free stream Mach Numbers is computed using 
the Karman-Guderley transonic approximation The location 
f the terminating shock wave as well as the pressure distribution 


s computed 


INTRODUCTION 


N THE PRESENT PAPER we Shall consider the flow over 
| a two-dimensional finite wedge at zero incidence for 
Mach Numbers in the lower transonic region (subsonic 
free-stream Mach numbers). 

In the analysis, the usual transonie small perturba 
tion method! is used. The resulting flow differential 
equation in the physical plane is still nonlinear, but 
what otherwise appears to be a well-defined boundary 
value problem arises if in the domain of the problem, 
There exist, 
however, Nikolski 


Taganov? as well as extensive experimental evidence 


the solution is assumed to be continuous. 
theoretical arguments by and 
which show that a continuous solution for the present 
problem is impossible, and that the solution must con- 
tain a prescribed discontinuity in the form of shock 
waves at an unknown location within the flow. 

In the analysis, the hodograph transformation will be 
made to linearize the flow differential equation; a 
wedge profile is selected in order to obtain, as far as 
possible, known boundaries in the hodograph plane. 
The difficulty due to the discontinuity mentioned above 
manifests itself in the hodograph representation in that 
a part of the boundary and the condition to be pre- 
scribed there are unknown, both being dependent upon 
the solution itself. In other words, when we formulate 
the problem in terms of explicitly defined mathematical 
problems we find that we have a transonic boundary 
value problem and a supersonic initial and boundary 
value problem. Part of the boundary condition in the 
former problem depends upon the solution of the super- 
sonic problem, while the initial data for the latter must 
be obtained from the transonic problem. 

Thus, in either the hodograph or the physical planes 
the presence of the shock wave prevents the straight- 
forward treatment of the problem. 

The first consideration to the present problem was 


given by Cole,* who gave a representation of the free- 
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stream singularity but avoided the difficulty described 
above by considering only a certain portion of the sub 
sonic boundary conditions. (The supersonic condi 
tions include the requirement that the wedge shoulder 
be on the zero streamline and the jump conditions con 
nected with the shock For Mach Numbers 


close to one, the resulting drag values obtained by Cole 


Wave. ) 


are reasonable when compared at Mach Number one 
with the results of Guderley.’ It is not clear, however, 
how the disregarded supersonic conditions affect the 
the drag and the pressure distribution along the wedge 
at still lower Mach Numbers. In_ the 


interior of the flow Cole's result leads to a physically 


free-stream 


impossible flow pattern. 

Trilling and Walker’ reconsidered the problem by 
fulfilling, in addition to the subsonic conditions, the 
condition in the supersonic region at the wedge shoulder. 
The flow pattern obtained by them appears to be un 
reasonable despite an allowance made for their dis 
regard of the shock conditions. 

In the present paper the complete solution is ob- 
tained by an iterative procedure in which all boundary 
and jump conditions are fulfilled. No rigorous con- 
clusion on the convergence of the procedure is at 


tempted, however. 


PHYSICAL DESCRIPTION OF THE FLOW FIELD 


Let us first consider the flow pattern over the wedge 
in the physical plane (Fig. la). Here the free stream 


of Mach Number 1/)* (based on the critical velocity) 
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and zero inclination relative to the wedge chord ap- 
proaches the wedge from the left. At the nose of the 
profile at point B, we have a stagnation point, while at 
the shoulder at point CD we have a Meyer expansion 
with the upstream condition being sonic. The expan- 
sion around the shoulder creates a supersonic region 
which is bounded on the upstream side by the sonic 
line CGF, while it is terminated on the downstream 
side by a shock wave FE’. The location of this shock 
wave is not known in advance. The flow external to 
this region is subsonic. 

In Fig. la the maximum ordinate of the supersonic 
region is shown at the juncture of the sonic line and 
the shock wave. This is contrary to the experimental 
results of Wood and Gooderum® and Bryson’? who 
found with an interferometer that the maximum of the 
supersonic region occurred along the sonic line itself 
but not at its juncture with the shock wave. The 
former configuration of the sonic line was, however, 
assumed in order to avoid streamlines along which the 
Mach Number first then decreases 
through one without a shock wave. 

The hodograph plane corresponding to Fig. la is 
next shown in Fig. 1b. In the hodograph plane the 
abscissa is essentially the difference of the Mach Num- 
ber from one, while the ordinate is the local streamline 
inclination with respect to the free-stream direction. 

In Fig. 1b point A represents the free-stream condi- 
tions. Since all streamlines must originate and _ter- 
minate at this point, it will be a singular point of the 


increases and 


stream function. 

The course of various streamlines is also illustrated 
in Fig. lb. It is seen that some of the streamlines de- 
part from and then return to point A without leaving 
the subsonic region; while others enter the supersonic 
region, encounter the shock wave at some point along 
the line FD’E, jump back to the subsonic region to a 
point along E’F, and then re-enter the free-stream sin- 
The location of the loci of the 
namely, FD’E 


guarity at point A. 
conditions before and after the shock 
and E’F is not known in advance. 
The shock wave at its juncture with the body at 
point E must be a normal shock. Thus, according to 
the transonic shock relations, points E and E’ in the 
hodograph plane will be equidistant from the sonic line. 
In Fig. 1b point E’ is shown at a higher Mach Num- 
ber or to the right of the free-stream singular point at 
That it cannot be otherwise can be seen 
assume that the relative 


point A. 
from the following reasoning: 
position of the two points is reversed such that point A 
lies closer to the sonic line than point E. Then let 
the free-stream Mach Number decrease toward zero 
that is, let point A move to the left. Since there is a 
minimum value of the Mach Number which point E’ 
may attain, there will then be one Mach Number for 
which points A and E’ would coincide, which is highly 
improbable. 

It is also to be mentioned that the region EDJD’E 
in the hodograph is on a second sheet other than that 
of the subsonic region, while in the physical plane the 
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flow is single-sheeted. This double-sheeted configura. 
tion of the flow pattern arises when a Meyer expansion 
occurs adjacent to a subsonic flow field. 


MATHEMATICAL FORMULATION OF THE 
BOUNDARY VALUE PROBLEM 


In the present paper we shall be concerned with the 
flow of a perfect nonviscous gas which fulfills the 
quirements necessary for the transonic approximations 
to be valid.'| With this approximation the basic floy 
equation for the dependent variable, the Legendre 


potential ¢, is given by 
Yun — NG = V | 


where the subscripts denote partial differentiations 
Here the independent variables y and @ are defined by 


n= A (2 Dy rv 1) Sy aad 1), 64, = 


where .\/* is the local Mach Number based on the criti- 
cal velocity a* at which the flow velocity is equal to 
the sound velocity, and 0 is the local flow inclination 
6) is the semi-nose angle of the wedge, and y is the usual 
ratio of specific heats. 

The transformation equations from the hodograph 
to the physical planes are given by 


(1 72) « 
> == al 
ll ae | 


x = a*(y+ 1) 

y = a*hoyv = f 

where ¢ and ¥ are the Cartesian coordinates in the 

physical plane with the Z-axis being parallel with the 

free-stream direction; x and y are reduced coordinates 
which we shall use later. 

Consider next the boundary conditions. Along the 
line of symmetry upstream of the wedge as well as 
alone the wedge surface——that is, in the hodograph 
plane along ABCD, along DE, and along E’A 
ge = 0. The location of points E and E’ is, as men 
tioned previously, not known in advance. At a far dis- 
tance from the body the flow approaches a_ parallel 
Mach Number J7/o* and zero flow 


we have 


stream with the 
inclination. 

By hypothesis, the free-stream singularity at point 
A is a logarithmic one as in the case of an incompressible 
flow over a semi-infinite body. 

Lastly, we must fulfill the jump conditions 


¢ge(m, A) = ¢e( Ne, G2) 


(4a 
¢y(m, 91) = ¥,(12, 2) § 


from a point (m, 4:) along FD’E to a point (m2, 42) along 
FE’, where 

(dg, dg6)su = (m — 2)/(02 — 0) | 
. (4b 
) 


6. — 0, = +(m — 2)V (m + 22)/2\ 


and the subscript sw identifies the bracketed quantity 
with the shock wave. 

Summarizing, our problem is now to find a solution 
of Eq. (1) in the domain ABCDED’FE’A which pos- 


sesses a logarithmic singularity at point A, whose pat- 
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tial derivative with respect to @ vanishes along ABCDE 
and E’A, and which fulfills the jump conditions of Eq. 
ja) and Eq. (4b) from ED’F to E’F. The solution 
must in addition be continuous and possess continuous 
first partial derivatives everywhere in the domain of 
(he problem with the exception of the free-stream singu- 


lar point. 


METHOD OF SOLUTION 


An analytical treatment of the boundary value 
problem given in the previous section is complicated, 
among other things, by the presence of the double- 
valued portion of the supersonic region. In order to 
avoid this difficulty and to obtain an explicitly defined 
problem of a standard form, we shall separate the con- 
siderations of the subsonic region from the supersonic 
and compute the latter region by the method of charac- 
For this purpose consider the hodograph 


Here if the location of the 


teristics. 
plane shown in Fig. Ib. 
boundary E’F and the condition to be prescribed there 
are known, then the solution in the entire subsonic re- 
gion (more exactly, the region ABCHFE’A) is com 
pletely determined. This then yield 
initial conditions along the sonic line which together 


solution will 


with the previously prescribed boundary and jump 
conditions will be sufficient to compute the remainder 
of the flow including the shock wave terminating the 
supersonic region. 

In order to find this unknown boundary condition 
we shall start from the solution obtained by Cole* which 
possesses the required singularity and which fulfills 
the conditions along the side of the wedge and along the 


axis of symmetry of the flow. This solution is given by 


I 9) /9)\1/35 — (1/6 1/2) pp , 1/3) /2)] 
gy’ = — (2/3)°'°2 T iT(1/6)/T(2/3)| X 

1/6 io 1/6 , J - ? o/s 2 ~ 

wu >. F(1/12, 7/12, 2/3; ¢, (5 


where F is the hypergeometric function, I is the gamma 


function, 
- 3/2 
w = (2/3)(—7n) 
On = 1(@ — 2n)* + w* T wy” | Zw 


and w) is the free-stream value of w. In Eq. (5) the 
coefficient of ¢''’ has been selected so as to obtain a 
unit length of the wedge nose. 

The logarithmic nature of the singularity at the free- 
stream position (o = 1) can be verified by an analytic 
continuation of the above hypergeometric function to 
the vicinity of ¢ = | by using the Barnes contour in- 
tegral representation of the hypergeometric function 
and noting that the difference of the resultant odd hy- 
pergeometric parameter from the sum of the two inter- 
changeable ones is zero. 

The condition gs = O along the characteristic CH 
is next fulfilled by a superposition of regular solutions 


g®? = » ® a,lT,(n) cos nn6 (6) 


n l 


where 
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H(n) = 3 \n)" "HH 3"(inrw) 1 < 0 

= 37 19/T(2/3) |(nx)~" 9 = 0 

= (2/3'/*)y' “[Jij3(nm@lw!) + J_1/3(nm!\w!) 


n> 0 


Here the usual notation for the Bessel functions is used 
and w = (2 3)(—7)’ “. The coefficients a, are de 
termined such that gy‘? + gy”? = 0 along CH by the 
method of least squares. 

In the physical plane the sonic line, as determined 
from the sum of the solutions ¢"'’ + ¢"” obtained above, 
is incompatible with the shock conditions. Thus, we 
must seek a regular solution ¢“, which will correct the 
sonic line by distorting the y-coordinate in its depend 
ence upon @ by an amount Ay(@) = ge’(0, @). In the 
selection of the Ay(@) it is first observed that the por- 
tion of the sonic line near the body as given by ¢''’ 4 
¢” should be essentially correct and unaffected by the 
nonfulfillment of the shock conditions. Secondly, it 
is to be noticed that locating point E in the supersonic 
region will locate the characteristic GJE, and conse 
quently point G on the sonic line. (Here it is to be 
remembered that points E and E’ are equidistant from 
the sonic line. Moreover, since point A is located near 
the sonic line, an approximate location of point E’ can 
be made without excessive error.) The correction 
function Ay must now be chosen such that in the con 
struction of the supersonic region CGFD’EDC by the 
method of characteristics the location of the Mach 
wave GJE is such that the »x-coordinate of point E 
is equal to that of point F. (Here the terminating 
shock wave in the first approximation is assumed to lie 
along a line of constant wx.) 

With the selection of a Ay(@), the supersonic region 
can be computed by the method of characteristics, and 
the shock wave can be constructed by an iterative 
method from point F using the shock conditions (4b) 
and assuming as a crude first approximation that the 
value of @, behind the shock—the subscript 2 denoting 
the conditions downstream of the shock 
linear function of y. (We shall see later that the exact 
choice of this function is unimportant, within limits, 
for the determination of the shock wave in the physical 
plane.) In this way the first approximation of the sub- 
sonic boundary E’F is determined. At this location, 
the x-coordinate from the solution formed by ¢"' 
¢’ can be computed; this will be in contradiction to 
the distribution of x(2, 62) from the method of charac- 
This difference in the x» will then be the con 
along E’F for the subsonic problem. 


is a suitable 


4 


teristics. 
dition for Ly *) 
Along the sonic line CGF we shall prescribe the change 
in the x-coordinate Ax(@) = ¢, (0, 0) due to the Ay(@) 
from the method of characteristics. This is equivalent 
to fulfilling the condition gy’ = 0 along the charac 
teristic of the Meyer expansion by a procedure of itera 
tion. 

With these boundary conditions and the two previous 
conditions along E’AB and BC, one can determine 
¢’ in the remainder of the subsonic region by the re 
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laxation method. Here, instead of the usual pointwise 
relaxation, we shall use a stripwise relaxation process 
in which we replace only the derivatives in one direction 
(the @ direction) by finite differences. In this way a 
system of simultaneous ordinary differential equations 
in 7 is obtained as an approximation to the partial dif- 
ferential equation. The further procedure is straight- 
forward; the reader is referred to reference S for de- 
tails of the procedure. 

The above subsonic boundary value problem can also 
be solved in principle by a series solution of the form 
given in Eq. (6); however, the number of terms re- 
quired for a satisfactory solution is so large that only 
an approximate solution is found, which will be used 
as the zero approximation for the relaxation method. 

Using the new distribution ge(0, 6) resulting from 
the above subsonic computations, ve now determine 
the second approximation to the sonic line and the 
supersonic region by the method of characteristics. 
The shock wave is now reconstructed using Eqs. (4a) 
and (4b), and, instead of assuming that the y-coordi- 
nate on the subsonic side of the shock is a linear function 
of @ as before, we locate the second approximation of 
the boundary E’F such that the y-coordinate as found 
from the method of characteristics agrees with the y- 
1) + yg” + 3 3) 

Again the change Ax along the new shock location 
E’F and the Ax along the sonic line due to the new 
Ay(@) will form the new conditions for the subsonic 


(3) 


calculations for the next higher approximation to yg” 


coordinate due to the solution ¢' 


In this way we now first fulfill the condition of the 
Meyer expansion at the wedge shoulder. This condi- 
tion can be fulfilled even though the conditions at the 
shock are not yet completely satisfied, since it can be 
expected that perturbations from the downstream side 
of the shock (except in the neighborhood of the juncture 
of the shock with the sonic line) will not influence the 
sonic line itself. 

Thus, at any point in the procedure we have all 
bourdary and shock conditions fulfilled except either 
of the two parts of the shock condition (4a)—that is, 
either the x(7., 0.) or v(m, 02), behind the shock obtained 
by the method of characteristics and that obtained 
are incompati- 


(3) 


directly using the subsonic solution ¢ 
However, from the boundary condition yg = 
we Call ex- 


ble. 
0 along the line @ = O—.e., 
pect that the x-coordinate near the y axis will not vary 
with @. Indeed, from preliminary computations of the 
subsonic region we find this behavior of the solution 
throughout the subsonic domain where the downstream 
conditions of the shock wave may be expected. More- 
over, from the supersonic computations it is also found 
that the position of the shock wave is practically inde- 
pendent of the value of @, within the expected range of 
@.— that is, the coordinates x and y of the shock wave 
are dependent explicitly on yn, (and of course on m) and 
the value of 7, as determined from the shock condition, 
except in the immediate neighborhood of the sonic line, 


Gyo = Xo = 0 


does not vary with the value of .. 
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As a consequence of this property of the solution 
after the condition of the Meyer expansion is fulfilled 
and the compatibility requirement on the x-coordinate 
is satisfied by the solution of the subsonic relaxation 
process, we can proceed to construct the shock waye 
by fulfilling the remaining requirement on the 
coordinate at the shock as previously described without 
essentially redisturbing the match of the x-coordinates. 

In other words, the shock wave can be reconstructed 
with the proper matching of the y-coordinates behind 
the shock, and the value of 7, for a given point on the 
shock in the physical plane will be essentially un- 
changed —-@, of course will change. Now, at the new 
location of the boundary E’F in the hodograph, the 
required boundary condition for ¢,'*’ will still be satis. 
fied because of the property that along E’F the de 
pendence of Av = ¢,'” 
it is practically independent of @. 
remembered that the remarks concerning the behavior 
in the neighborhood of E’F also 


on 7» is almost unchanged, and 

(Here it is to be 
(2) (3) 

Ty + 2 


', and the contribution of the latter to the 


of ¢' P 
: (3 

pertain to ¢ 

x-coordinate predominates over that of the other two. 


_ ‘ ‘ ° 3 p — 
rhus, the prior approximation of ¢” will fulfill ap 





proximately the required condition, and a satisfactory 
value of ¢ can be expected after repetition of the 
above procedure. 

The solution obtained in this way is now to be com 
and ¢'”’, and the required results in the 
physical plane are then obtained by using the transfor- 


mation (3). 


DETAILS OF THE COMPUTATIONS AND 
DISCUSSION OF THE RESULTS 


The method described in the previous section has 
been applied to compute the flow pattern for the value 
of the (2/3) " = 0.5. 
Here the transonic parameter is defined by & = (1 - 
Mo*?) [Cy + 100/77? \Jo* is the free stream 
Mach Number based on the critical velocity, is the 
half angle of the wedge, and y is the usual ratio of the 
In order to obtain the influence of the 


transonic parameter a) = 


where 


specific heats. 
transonic parameter on the drag coefficient, the surface 
pressure distribution along the wedge was computed 
for the additional value of w = 0.5. Here for the 
computation of the drag (as we shall see shortly) only 
the portions of the solution ¢'’’ and ¢'”’ are required 

For the calculation of ¢‘~’ required to fulfill the con- 
dition at the wedge shoulder, ten terms of the series 
given in Eq. (6) were found to be sufficient (see Fig. 
2). For the computation of the relaxation portion of 
the solution, the interval between the strips was taken 
successively as A@ = 0.05, and 0.025, and in critical 


regions through the shock wave; also the smaller inter- 





val 0.01 was used. 

The resulting values of ¢’ are shown in Fig. 3. 
It is to be noticed here that ¢'”’ attenuates very rapidly 
as the surface of the wedge is approached at @ = |, 
having negligible influence at that location. Thus, 
for the computation of the pressure distribution along 


(3) 
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Construction of the supersonic region and a comparison 
of the sonic lines from various approximations 
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the sloping side of the wedge, and, therefore, of the drag, 


only the knowledge of ¢'’ and ¢™? is required. 

The resultant calculation of the supersonic region by 
the method of characteristics in next shown in Fig. 4, 
giving the location of the terminating shock wave and 
the conditions before and after the shock. Here a 
solution was obtained with the maximum ordinate of 
the supersonic region occurring at the junction of the 
sonic line with the shock wave. However, there was 
no evidence in the numerical calculations that definitely 
excluded the possibility of other configurations. Also 
shown in Fig. 4 is a comparison of the sonic line from 
the various approximations. Here one clearly sees 
the inadequacy of the various approximations in giving 
the correct location of the sonic line. 

In order to compute the pressure distribution over 
the wedge let us recall the transonic form of Bernoulli's 


equation, 
Cp* = (p — p*)/q* = —2Ay + 17 O"y (7) 


where C,* is the pressure coefficient relative to the 
sonic conditions, p is the static pressure, g is the incom- 
pressible ‘“dynamic’’ pressure, and the asterisk denotes 
the sonic value of the latter two quantities. In Fig. 5 
we give the resultant pressure distributions together 
with that for Mach Number one from referetice 4. 
Also shown here is a comparison of the pressure distri- 
bution resulting from Cole’s approximation with the 
present results for w) = 0.5. The agreement between 
the two is good and is sufficient for engineering purposes. 

The drag coefficient can next be found by integrating 
the streamwise component of the pressure coefficients 
along the wedge surface. In Fig. 6 we show a compari- 
son of the present results with the experimental results 
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of Bryson’ as given in Fig. 7 of Spreiter.’ Here the 
experimental data have been plotted according to the 
parameters suggested by Spreiter®; these parameters 
differ from those resulting from the usual transoni 
theory of Guderley and von Karman by higher order 
quantities, but their use results in a considerably better 
correlation between theory and experiment. It is seen 
in Fig. 6 that the present results lead to an improved 
agreement with experiment as compared to Cole's re 
sult. 
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Loca range of slope ratios, leading-edge slope to Mach wave slope, 
Flow, from 0.15 to 4.1. Lift or normal-force characteristics for these 
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eed presented as normal force coefficients per radian for the Mach 
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‘edg with the linear theory both as a function of Mach Number and 
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DESIGN RANGE 
RANGE OF TESTS COMPLETED 


20 25 3.0 3.5 40 
tonw = tonw/tony 


Fic. 2a. Range of tan w/tan uw for indicated delta wing models. 


INTRODUCTION 


b ive INVESTIGATION was undertaken principally to 
obtain continuous measurements of lift or normal 
force for a family of delta wings (see Figs. la, |b, and Ic) 
over the Mach Number and Reynolds Number ranges 
that could be obtained in the 12-in.' and 20-in. super- 
sonic wind tunnels* at the Jet Propulsion Laboratory. 
The test Mach Numbers were selected such that force 
measurements for any one wing could be checked against 
another for a range of the parameter tan w tan u 


BAR 4. The Mach Numbers plotted against tan 
w/tan p BAR 4 is indicated in Figs. 2a and 2b. The 
440 
4.00 + + + + + 
3.60} + } | + + 
320 } | | | | | | 
2.80) + + + + + + + 


tonw/ tony 














effects of Reynolds Number on the force characteristics 
of these wings were studied in the 20-in. wind tunnel 
A somewhat similar program was carried out by Love 
for similar wing profiles over a much more restricted 
range of Mach Numbers. 

Force measurements were made for increments of a 
0.2” in the range —1° <a < 1 For angles of attack 
outside this range, @ was increased in steps of | and 2 
degrees to check the linearity of the theory at the higher 
angles.” For certain wings, the angle of attack was al 
lowed to vary up to the largest possible values of a for 
which the wing components could safely be expected to 


withstand the bending loads. 
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MACH NUMBER 


FiG. 2b 


Slope ratio tan w/tan w vs. Mach Number for delta wings tested in 20-in. wind tunnel. 
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NORMAL-FORCE CHARACTERISTICS OF DELTA WINGS 669 
THEORY = 
PABLE | 
. . , : Geometric Characteristics of Delta Wings Teste 
Normal-force measurements were made as a function silti eee elta Wings Tested 
of angle of attack* for all wings listed in Table | at the ~~ > ss 
Mach Numbers and Reynolds Numbers indicated in | és + o86 
a Gea | A 
Figs. 2 and 3. To facilitate the comparison of these ie : n) : 
measurements with linear theory*® it is convenient to ‘ 
I ‘ vs : : d . or t => -/\ 
express these results in the form of normal-force co I \| 
os se .. | 1 
eficient per rad. Cy. a So hy Vy 
a = ie P : AS > + I, i.) 
In Figs. 4 and 5 the lift or normal-force slope ratio \ ma | 
. . *¢ e | Ht 
3 4)Cy. (4/8 theoretical-lift-curve slope for a delta \ % t, -f jt 
a \ 70 | /|' jt 
a . ° > ee 2 ; /\\| | 
wing of infinite extent) is plotted against the slope ratio / | \ , i/ ; \i 
| 
tanw tan yu (ratio of leading edge slope to Mach wave f- - 4 60 43741 y 
slope) for the families of blunt and sharp leading-edge * ee ay 
wings, respectively. In Fig. 6 the slope ratio (8 4)Cy 
- a 
is plotted against tan » tan uw for A,5, A.5, and A,» to ; pss p= ‘ueen Gee Pen EN - 
- “¢ ee: Configu- by | ' 
illustrate the effects on the lift due to the variation ot peace 6 * | * ~ i) | Ye , . 
he wi rofile for wines of oe , ‘ r 4 - t —— ——}————$_+—_—__— +——_—++- 
the wing profile for wings of the same plan form. It is ss | aan | 2450 | osm | ses2 | ass | «a0 | 1.069 | 0900 
sometimes convenient in the case of delta wings to com A,5 4.462 | 2.467 | 0.0575 | 3.616 | 0.904 | 42.11 | 1.063 | 0.950 
¥ . : AS 4.468 | 2.449 | 0.1225 3.632 | 0.908 42.25 1.225 0.956 
pute the normal-force coefficients using the span , | ome tate | dee | ame fome | wey | tan ue 
squared as an area parameter rather than the actual 4.3 | 2.932 | 4.093 | 0.0954 | 1.433 | 0.358 | 19,70 | 1.967 | 1.520 
ns are: % thi 4 ts c ] : tl : rasr ot . \,3 2.932 4.093 0.0954 1,433 0.358 19.70 1.967 1,520 
wing area. oelicients computer 1 Hs way five a 42 2.275 5.304 0.1232 0.858 0.214 12.08 2.765 1.985 
somewhat clearer picture of the variation of lift for \,2 | 2.272 | 5.298 | 0.1232 | 0.858 | 0.214 | 12,08 | 2.765 | 1.985 
P - . Al 1.769 6.831 0.1592 | 0.518 | 0.129 | 7.35 4.035 | 2,600 
values of the slope ratio tan » tan w approaching zero. : jen | | | ‘& 
Fig. 7 illustrates the case where the lift or normal-force 
slopes are computed in this fashion and are plotted 
against the slope ratio for all wings tested. In Figs. 8 ing-edge case on the basis of the slope ratio namely, 


and 9 the lift or normal-force slopes are compared with 
theory as a function of Mach Number for the various 
wings tested. 

In the discussion which follows the wings will be 


classified either as the subsonic or the supersonic lead- 


whether tan w, tan yu is less than or greater than unity. 
Although the experimental measurements which will 
be discussed indicate no sharp demarcation between 
these regimes, the linear theory with which these results 
are compared makes such a distinction and hence pro- 
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vides a way of talking about the experimental results from the wing apex are far removed from the leading 


categorically. 
RESULTS AND DISCUSSION 


Subsonic Leading Edges 


This is the case where the slopes of the Mach waves 
in the flow are normally greater than the slopes of the 
For this 


leading edges of the delta wings in question. 
case the theory predicts that the edges of the wing will 


experience a subsonic-type compressible potential flow 


where the component of the free-stream velocity normal 


to the edge is subsonic. It may be noted in Figs. 4 and 





wing conform 
the 
appear in 


edge namely, where conditions at the 
closely with the assumptions of the linear theory 
experimentally measured values for (3/4)Cy_, 
general to be slightly greater than but in substantial 
agreement with the theory, that is, for the range 0.14 < 
tan w/tan uw < 0.60. Good agreement also may be ob- 
served in the pressure measurements which were made 
by Drougge.® For values of the slope ratio greater than 
0.60-—-that is, where the linear theory curve appears to 
intersect the data, through tan w/tan pu 1.0 the ex- 
perimental values fall between 10 and 15 per cent below 
the linear theory for wings with blunt and sharp leading 
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the subsonic and supersonic regimes of 


ground tan w 
tional point insofar as the experimental normal-force 


measurements 


NORMAL-FORCE CHARAC 


elect of the leading edge, whether sharp or blunt, ap- 
neared to have little effect on the variation of (8 Cy, 
the 
| appear to sustain a maxt- 


with tan w/tan yw. However, wings with blunt 
leading edges shown in Fig. 


num of approximately 3 per cent more lift than the 

sharp-edged wings over the same range of slope ratios 

1)Cy 
a 


chown in Fig. 5. The maximum deviation of (& 


above the theory for this range of slope ratios 0.14 < 
tanw tan w < lis less than 10 per cent. Significantly, 
the region in the vicinity of the junction point between 
the theory 

tan pu | does not appear as an excep- 


ire concerned. It should be pointed out 
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that, while the Mach waves of the flow have the same 


or nearly the same slope as the leading edge, the shock 


waves from the wing apexes are still appreciably ahead 


of these edges, and there may yet be substantial interac- 


tion between the upper and lower surfaces of the wing 


around the leading edges (see Figs. 10 


and 


11). 


The 


wings in fact are beginning to experience flows, locally 


at least, which could be considered similar to those that 


would be expected along the edge of a two-dimensional 


wing in transonic flow.” 


s 


The agreement of the experimental measurements 


with theory as a function of Mach Number may be seen 


in Figs. S and 9. 
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MACH NUMBER 


Supersonic Leading Edge 


In this range, where the slope of the Mach waves of 
the flow exceeds the slopes of the delta wing leading 
edges, the linear theory presupposes that local wave at- 
tachment has occurred all along the wing leading edges 
and that the component of the free-stream velocity nor- 
10 and 


11 it may be seen that, for the wings tested here, 


mal to the edge is supersonic. However, in Figs. 
there 
exists a substantial part of the range tan w tan u > | 


for which the shock wave from the wing apex still ap- 







Normal-force coefficient per radian vs. 





Mach Number for wings indicated 

pears to be detached from the leading edges. Further 
more the shape of the wing edges, whether sharp or 
has considerable effect on the local flow condi- 
Specifically, the A,5 and A,4 schlieren photo- 
10 and 
attachment did not occur throughout the range of test 


blunt, 
tions. 
that the local wave 


graphs (Figs. 11) indicate 


Mach Numbers. The flow locally around such edges, as 


indicated previously, might be classified as transonic, 
and qualitatively at least one would expect the flow 
conditions normal to these edges to display characteris- 


tics similar to those experienced along the edges of two- 
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NORMAL-FORCE 


dimensional airfoils with detached bow waves. A large 
portion of the lift of a delta wing is carried or sustained 
along its leading edges; hence, the edge effects tend to 
dominate the picture insofar as lift measurements are 
concerned. In Fig. 10, for instance, as the shock wave 
from the A,5 wing apex approaches the edge of the wing 
the departures from theory are particularly apparent; 
the normal-force slope ratio (@ 1)Cy, in Fig. 4 falls off 
sharply in the range | < tan w tan w < 1.5, then re- 
covers up through tan w tan uw & 2 where (8 Cy, is 
about 90 per cent of linear theory. At this value of the 
slope ratio, as may be observed in Fig. 10, the plan form 
schlieren photographs indicate that the position of the 
shock wave is almost parallel to the leading edges and 
appears to retain its same relative position to the edge 
throughout the test range.. The variation in (8 HCy, 
with slope ratio tan » tan w for values larger than 2 
appears to increase almost linearly up to the largest 
namely, +.04S8 where 
Whereas the 


value of the slope ratio obtainable 
84)Cy, is 97 per cent of the theory. 
4,5 wing is notably affected by the motion of the de- 
tached bow wave relative to its edges, the A,4 wing lift 
appears to be more orderly and is affected to a lesser de- 
gree (see Fig. 4). This wing sustains more lift-—5 to 15 
per cent over an identical range | < tan w/tan uw < 
275 than the 4,5 wing. The agreement with theory 
for 4,4 improves greatly with Mach Number, as can be 
observed in Fig. 8, where at the highest test Mach Num- 
ber the measured normal-force slopes are 2 to 3 per cent 
below the theory. It is apparent that the differences 
exhibited by these two wings as a function of tan w tan 
u cannot be resolved by a simple force measurement 


test; the effect of edge bluntness and wave attachment 
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Fic. 10. Position of shock 


CHARACTERISTICS 


1.97 
1.55 











OF DELTA WINGS 6 


requires a more extensive investigation which should in- 
clude a study of viscous effects. However, as will be 
demonstrated, certain edge effects have been obtained 
by studying three wings of identical plan form but with 
The wing singled out in that in 
Before 


different wing profiles. 
vestigation was the Ad wing (see Fig. 6). 
elaborating on this phase of the investigation, the be- 
havior of the sharp leading-edge delta wings in the range 
tan w tan uw > | will be discussed. The wings designated 
as A, are the sharp-leading-edge versions of the blunt- 
leading-edge counterparts. The leading edge of these 
wings is a 12° wedge in the stream direction. These 
wings are also affected by the close proximity of the bow 
wave to their edges, but the effect is not so pronounced 
as for blunt-edged wings. In the plan-form schlieren 
photographs of Fig. 10, wave attachment does occur for 
the A,5 wing but at a much higher Mach Number than 
assumed by linear theory (see Figs. 5 and 6). Before 
attachment occurs, however, it may be observed in Fig. 
5 that the A,5 wing in particular sustains from 10 to 15 
per cent less lift than predicted by the linear theory. 
Comparing the plan-form schlieren photographs in Fig. 
10 with the normal-force characteristics of Fig. 5, it 
may be noted that after wave attachment the lift pre- 
dicted by theory is exceeded by about 6 per cent for the 
AS wing. A simple method employed for predicting 
wave attachment which was in agreement with the ob- 
served results of the plan-form schlieren photographs 
was to determine the minimum value of tan w tan u for 
wave attachment to a two-dimensional wedge’ having 
the same wedge geometry and Mach Number as the 
delta wing in a direction normal to the edge. The re- 
sults of these computations are shown in Figs. 5 and 6. 
In Fig. 5 the experimental results for the 4,5 wing de- 
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MACH NUMBER 1.79 1.97 2.15 2.34 2.55 2.79 3.05 3.35 3.67 4.06 
tanw/tanz=B asa 88 1.01 1.14 1.28 1.41 1.56 1.73 1.85 2.12 2.36 
ELLIPTIC LEADING EDGE 
A.4 
RADIUS / 
w= 30.9° 
SECTION A-A 
tanw/tanu 7 ; : ; - 1 
<a 
MACH ANGLE | MACH ANGLE | MACH ANGLE | MACH ANGLE | MACH ANGLE | MACH ANGLE | MACH ANGLE | MACH ANGLE | MACH ANGLE | MACH ANGLE| 
Hi 349° | 305° | prar7? | s283° | sear? | pe2io® | psio2® | yst79® | prise? | pcre ze | 
Fic. 11. Position of shock waves relative to leading edge of A,4 wing. 


parts noticeably from those of A,5 at tan w/tan wp = 1.2; 
the A,3 wing appears to sustain about 8 per cent more 
lift than A,5 at the largest values of the slope ratio ob- 
1.59. Although 


the A.3 wing's leading edge was a 12° wedge in the 
: § 8 


tainable for this wing, tan w tan p = 


stream direction, the section in the direction normal to 
the edge was 34°, which is twice the wedge angle of the 
A.5 wing. The normal-force measurements as a func- 
tion of the slope ratio for A,3 appeared to be in better 
agreement with the blunt A,4 wing than the A,5 wing 
as may be observed in Fig. 8. In Fig. 8 where Cy, is 
plotted as a function of Mach Number, the lift sustained 
by the A,3 wing, similar to the 4,4 wing, appears to be 
in better agreement with the theory at the highest 
Mach Numbers where the indicated measurements are 
5 to 5 per cent below the theory. 

The notable differences in the behavior for Cy, be- 
tween wings over the supersonic-leading-edge range 
may be ascribed to many factors which are not directly 
evaluated in this investigation, such as the effects of 
viscosity on flow separation along the edges as a function 
of wing thickness and edge bluntness. However, cer- 
tain general results and conclusions may be reached in- 
sofar as the geometry of the leading edge is concerned 
by comparing the sharp- and blunt-leading-edge wings 
of the same plan form. In Fig. 6 such a comparison has 
been made for wings of identical plan form, over the 
range of tan w/tan u > |. In addition to the A,5 and 
A.5 wings already discussed, a third wing, A,5, is in- 
cluded which had a symmetrical double-wedge profile 
in the stream direction (see Table 1 for details of its 
geometry). Furthermore the shock-wave configurations 
for the three wings in question have been compared in 
Fig. 10 where the A,5 wing, just as the A,5 and A,5 
wings discussed previously, is affected by the close 
proximity of the bow wave to its edges but to a lesser 
degree than either of the two flat-sided wings. It may 
be further noted in Fig. 6 that the wave attaches to the 
edges of the wing at very nearly the value of tan w tan 
u calculated by wedge theory. Prior to attachment this 
wing sustains from 5 to 10 per cent less lift than theory, 
whereas after attachment occurs the theoretical lift is 


exceeded, the excess increasing with increasing Mach 
At the highest test Mach Number, the lift 
The same general 


Number. 
was 10 per cent over theoretical. 
effect is observed in the case of the A.5 wing. As a rule 
of thumb it might be said for delta wings of the same 
plan form that, prior to wave attachment to the leading 
edges, the wing in the range tan w tan u > | sustains 
less lift than linear theory; whereas, after attachment, 
the theory is exceeded to a greater or lesser extent de- 
pending on the wedge angle of the edge. If the wave 
does not attach, linear theory will always overestimate 
the lift. 
or underestimating the lift, indicate the trends ex- 


All results, however, aside from overestimating 


hibited by the linear theory at the largest values of 
Mach Number and tan w tan uw measured, as may be 
observed in Figs. + through 9. It should be pointed out 
here that force measurements for both the A,5 and AJ 
wings were made at two different Reynolds Numbers 
(350,000 and 150,000 per in.) throughout the Mach 
Number range. It may be observed in Fig. 6 that for 
these particular tests there appeared to be no significant 
effect on Cy, due to Reynolds Number over the range 


indicated. There are however, comparatively large 
effects indicated by the variation in edge configuration 
It appears that for large values of tan w tan yu, say in 
excess of 2, the sharper the leading edge is made the 
more effective the wing becomes in producing lift. 
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SUMMARY 


The problem of the interaction between a plane entropy wave 


extent is investigated. Asa result of the interaction, it is found 


kinds of disturbances are present—namely, an entropy mode, 
, vorticity mode, and a sound mode. The nature of the sound 
wave generated depends on the orientation of the upstream dis 
Within 
sound 


turbance certain orientations of the upstream dis 


turbance, the Waves generated downstream attenuate 


Bevond these orientations, the sound waves generated have 
constant amplitudes 
not attenuated, there is no phase shift in the entropy disturbance 
across the shock; when the downstream sound waves are attenu 
shock. An illustrative example is given for the interaction be 
tween a plane normal shock and a sinusoidal entropy wave. At 
i given shock strength the amplitudes of the shock displacement 
ind the downstream disturbances generated are plotted as func- 


tions of the orientation of the upstream disturbance 


SYMBOLS 


U;, l velocity of the unperturbed main flow upstream 
and downstream of the shock, respectively 

A,, A speed of sound of the unperturbed main flow 
upstream and downstream of the shock, respec 
tively 

VW, M = Mach Number of the unperturbed main flow 
upstream and downstream of the shock, re 
spectively 

€ = inclination of the unperturbed shock plane with 
respect to the main-flow direction upstream of 
the shock, Fig. 1 

B = inclination of the unperturbed shock plane with 
respect to the main-flow direction downstream 
of the shock, Fig l 

p = pressure of the unperturbed main flow 

py = density of the unperturbed main flow 

C, = specific heat of the medium at constant pressure 

ratio of specific heats of the medium 

Ap, Ap . . 

As, Az = perturbation of pressure, density, entropy, and 

. velocity of the flow field 

= dimensionless quantities of the above perturba- 


tions as defined in Section (2), Eq. (2.1) 
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f arbitrary orientation and a plane oblique shock of infinite 


that, in the perturbed flow field downstream of the shock, three 


When the downstream sound waves are 


ited, a phase shift occurs in the entropy disturbance across the 
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Interaction of a Plane Shock and Oblique 
Plane Disturbances With Special Reference 
to Entropy Waves' 

CHE-TYAN CHANG* 
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coordinate system defined in Section (2) and 
Fig. 1 

time variable 

a reduced space variable corresponding to ¢, 
defined as r = At 

irrotational and rotational part of the perturbed 
velocity field 7, respectively 

scalar and vector potentials associated with 7 


, Eq. (2.3) 


and u 
perturbed vorticity associated with #, 
components of the vorticity perturbation # along 

the three coordinate axes ox, oy, and oz, respec 

tively 

shock displacement 

local shock deflection 

local shock oscillating velocity 

dimensionless perturbed velocity components 
downstream of the shock taken along ox* and 
oy*, respectively 

Mach Number corresponding to an equivalent 
normal shock upstream and downstream of 

the shock, respectively, Eq. (3.4 
shock strength as defined in terms of the ratio of 

the pressures of the unperturbed flow across 

the shock, Eq. (3.5) 
amplitude and wave number of the upstream 

entropy wave, Eq. (4.1) 
inclination of the upstream entropy wave with 

respect to the main-flow direction l’, Fig. 3 
directional cosines of the normal to the wave 


front of the entropy wave with 


respect to the coordinate system x,0y, 
drift speed of the upstream entropy wave along 


upstream 


the shock, Eq. (4.2), Fig. 3 


apparent main-flow direction downstream of 
the shock with respect to a mov ing observer, 
Eq. (4.3), Fig. 3 

inclination of U’, with respect to the downstream 
also inclination of 


main-flow direction Ul’; 


shear waves generated downstream, Fig. 3 
effective Mach Number and Mach angle cor 

responding to U’,, Eqs. (4.4) and (4.5) 
coordinate axes moving at a speed of C, along 


the shock, Fig. 3 


dimensionless perturbed velocity components 
downstream of the shock taken along OX and 
OY, respectively 
critical inclination of the upstream entropy wave 
corresponding to the drift speed (C,);, and 
C,)2, respectively, Fig. 4 
a dimensionless parameter related to the down 


defined as 


stream pressure perturbation, 
q=- p COS py 
distance measured along the shock with respect 


to the moving coordinate system y’ y* - 


(C,/A)r 
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Ay, a., 6, = amplitude of the downstream entropy wave, 
Fig. 7 
A,, dp, 6, = amplitude of the downstream sound wave at the 


shock, Fig. 8 


Ary, ag, be == amplitude of the downstream vorticity wave, 
Fig.9 

Ay, ay, by = amplitude of the shock displacement, Fig. 10 

d = decay distance associated with attenuated sound 
waves 

Ly = wavelength of upstream entropy wave 

rE. = projection of L; along the shock 


(1) INTRODUCTION 


Ws of previous investigators have shown that, 
for weak disturbances occurring in a uniform 
compressible flow, a decomposition into three modes is 
possible—-the vorticity mode that may be interpreted 
as turbulence in an incompressible flow, the pressure 
mode that represents a field of sound waves, and finally 
the entropy mode that may be viewed as temperature 
spottiness.! If the disturbances are assumed to be 
small with respect to the basic flow within the validity of 
first-order theory, these three modes are governed by 
independent differential equations. Interaction be- 
tween the various modes, however, can take place if 
there is a discontinuity surface, such as a shock, pres- 
ent in the flow field. 

As a result of the interaction with the shock, usually 
all three modes will be produced downstream for an 
upstream disturbance of a given mode. 

When the upstream disturbance happens to be prop- 
erly orientated with respect to the shock, it is possible 
for the shock configuration to be steady in time.‘ In 
most ‘circumstances, for disturbance of an arbitrary 
orientation, the shock configuration varies in time. 

So far as the upstream disturbance is concerned, the 
three modes may occur either separately or simul- 
taneously. Since both the boundary conditions and 
the governing equations of the perturbed flow field are 
linear, the principle of superposition holds; the in- 
fluence of the various modes occurring upstream on the 
flow field downstream, therefore, can be studied sep- 
arately. Shock interaction with upstream disturbance 
of either the vorticity mode or the pressure mode has 
been treated by Ribner? and Moore® in two separate 
NACA reports. 
stream disturbances of all three modes has been given 


A unified treatment concerning up- 


by the author in a preliminary report of Project Squid.* 
Since recently there is a renewed interest regarding 
shock interaction with entropy mode,'® an article deal 
ing exclusively with that subject may still be of some 
value for future investigations. 
Our analytical model now follows: A wedge is placed 
in the downstream side of a uniform flow field and an 
oblique shock is formed at the wedge when the up- 
stream flow strikes the wedge. The flow field thus is 
divided into two subregions an upstream region hav- 
ing a uniform velocity Ul, and a downstream region 
We now 


having a uniform velocity LU’ (see Fig. 1). 


imagine that plane waves of entropy disturbance are 
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introduced upstream and are drifting down with th, 
main flow toward the shock. 

Since our main interest in this paper is to study the 
interaction of the shock with the upstream disturbance 
and its effect on the downstream flow field, we shall 
purposely ignore the presence of the wedge from noy 
on. The shock is thus taken as infinitely extended 
so that reflection phenomena of sound waves generated 
downstream are neglected. (Readers interested jy 
these phenomena may consult the work of reference § 
As a further simplification, the medium is to be taken 


as a nonviscous ideal gas. 


(2) GOVERNING EQUATIONS 


For convenience of analysis, the flow parameters 
will be replaced by their corresponding nondimension 
Thus if Ap, Ap, As, and Ay denote th 
perturbations of pressure, density, entropy, and v« 


alized ones. 
locity, their corresponding dimensionless parameters 
will be denoted by 


Ap/ pm; 
s = As/C,, u = Anv/A (2.1 


p = Ap ypm, p= 


In the above, p,,, p,, are the pressure and density ot 
the unperturbed main flow, C, is the specific heat at 
constant pressure, y is the ratio of the specific heats of 
the gas, and .! is the speed of sound of the medium 
We shall also adopt the convention that, whenever 
no subscript is attached, the flow parameters are those 
in the region downstream of the shock; for similar flow 
parameters upstream of the shock, a subscript | will be 
used —for example, the pressure of the main flow up 
stream of the shock will be denoted by /,,,, while the 
pressure of the main flow downstream of the shock 
will be denoted by ?,,. 


Three sets of rectangular axes will be used: * and 


OX” ale 
oy*, with oy* taken along the shock plane; ox; and 04; 
with ox, taken along the velocity vector l; of the up 
stream main flow, and finally ox and oy, with ox taken 
along the velocity vector l’ of the downstream main 
flow. The relationship among these coordinate axes is 
indicated in Fig. 1. 

The equations governing the flow fields, both up 
stream and downstream of the shock, are the three con 
servation laws of mass, momentum, and energy. The 
equation of state, in addition, gives a relation between 
the three thermodynamic variables. 

If one decomposes the velocity field into two parts 
an irrotational part 1, and a rotational part 2, or, by 


definition, 


curl x7; = 0, div 4, = 0 (22 


one may then introduce two parameters—a scalar po 


tential ¢ and a vector potential defined by 
ui = —grad¢, us = curlé (2.0 


With respect to the two potentials, the governing 
equations can be written as the following: 
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Ad — (D*¢/Dr*) = 0 (2.4) 
divE =0, DE Dr =0 (2.5) 
DS/Dr = 0 (2.6) 


where the operator 

D/Dr = (0/07) + A,(0/0%) (2.7a 
for the flow field upstream of the shock, and 

D/Dr = (0/0r) + M(0/dx) (2.7b 
We have 


also replaced the independent time variable ¢ by two 


for the flow field downstream of the shock. 


reduced space variables, thus 
r= Al, T} = Ad = (4,/A)r (2.8) 


The three modes, entropy, vorticity, and sound, are 
clearly indicated by these equations. According to 
our nondimensional representation, the vorticity @ is 


given by 


® = A curl 4, = A curl curl £ (2.9) 
where © = 1é ‘is in 7 ke. Using the above condition 
of div & = O, this can be written as 

AAK = -G (2.10) 
S44 Ay, Ay Aig O Ss 
P+ A 1 A » Ags 0 p 
“,* Ag; Age Ass 0 nS 
v,* 00 60 Aalin* 


In the above, subscripts + refer to flow parameters immediately behind the shock; similarly, subscripts 
The downstream perturbed velocity Az has been resolved into 


flow parameters immediately ahead of the shock. 


The sound field is represented by the scalar potential 
¢ for 
p -—— — (Do Dr) 


ui = grad @ 


(3) BOUNDARY CONDITIONS 


By taking the mean position of the shock along the 


oy*-axis, the shock configuration can be expressed as 
x* = f(y", 1) (3.1) 


Within the accuracy of first-order approximation, the 
local perturbed velocity and deflection of the shock are 
given by .ly, (=y,) and y,«, respectively. 
tion of .l¥, is normal to the local shock front (see 
If we isolate a small element from the shock 


The direc 


Fig. 2). 
and superimpose a velocity vector of the same magni 
tude but of the opposite direction as 1y, to the whole 
velocity field ahead and behind the shock and then 
apply the usual Rankine-Hugoniot’s equations to the 
flow parameters across the shock, the downstream per 
turbed flow parameters can be solved explicitly in terms 
of the local shock property (its deflection and velocity) 
and the given upstream perturbed flow parameters;* 


thus 
Tit M cos 2 Tit 
wgirvesndel Oe” bed (3.2) 
tx, \J cos B 31 
4 lo 


reler to 


components Au* and Av* normal to and tangential to the unperturbed shock plane, respectively; u* and v* are simply 


their corresponding dimensionless forms, thus u* = Au* 
of the perturbed velocity upstream of the shock. 


the following: 


A and v* = Av*/A. 


The coefficients occurring in the operating matrices are given as 
5 5 


Similar notation applies to components 


An [(m/ Pim) (N/N,) |? — (y — 1) {1 (pm/ Pim) |N? 
Ay [(N2/(1 — N?)| (1 (Pm/Pim)| {1 + (4 1)N2] + $1 — [Com pim)?(N M1) ]?1) 
As) [VN (1 — N®)] [1 (pmV/pimNi)?| + [1 Pim) ly. N24 
Ay (4 1) [1 — (m/Prm)] [1 — (1/N12) (pm/prm) |N? 
\ —|[N?2 (1 N*)|] (Cl — (pm/pim) | + U1 (1/N12) (pm pim) | $1 + 4 1) {1 (Pm Pim) |N2!) 
\ [V/(1 N2)] 51 (pmN/ pimNi)?| + yf1 Pim) | [1 (1) N1?) (pm/ Pim) |N*} 
a vy — 1) [1 — (pm/Pim) ]2(N2/Ny as 
i [NV (1 — N2)] [1 — (pm /pim) | }2 + (y — DD — (pm! pim)|N2} (N/N, 
\ 1/11 — N2)] (N/M) 31 — [(pm/pim)N |? y{1 (pm / Pim) |? N24 
\ PmlV / pimi 
711 (y= Fe (Pim Pm) |? (Pm Pim) N 
Tr) -(N (1 — N®)] [1 (Pim Pm)|12 + (4 (pm Pim) |N?4 
131 (1 (1 N2)] [1 — (pim pm)]}1 + N2 + (4 (1 (pm / Pim) | N74 
141 [(pm/ Pim) — 1]N 
In the above, .V; and N are defined as 
N, = M, sine, N= MsinB (3.4 
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PLANE SHOCK AND 
and can be viewed as the Mach Numbers upstream and 
downstream of an “equivalent normal shock.” 

An examination of these coefficients shows that their 
yalues depend on the parameters Ni, N, and pp’ pim. 
If one introduces a parameter x defining the shock 


strength as 


| Pm Pim 


and takes the medium concerned to be air with y = 
1.40, it can be shown that 


N,? = (6x + 1)/7, N? = (x + 8) 7xI 


(3.6) 
Pm Pim = (6x + 1)/(x + 6) { 


This reveals that all the above coefficients are merely 
functions of the shock strength x. The obliqueness 
of the shock, or the dependence on the shock angle £8, 
enters only in terms involving Y,*«—1.e. the local shock 


inclination. Further examination shows that it af- 





fects only the first three equations of the system; this is 
due to the fact that, under the present notation, the 
term (.1/ cos 8)~,« simply means the contribution of the 
main flow to the oscillation of the shock due to its 


obliqueness (see Fig. 2). 


$) DESCRIPTION OF THE (GENERAL NATURE OF THE 


DOWNSTREAM FLOW FIELD AND THE METHOD OF 








ANALYSIS USED 


Since a plane entropy wave is characterized by its 
amplitude &,, its wave number &, and its inclination 6 
with respect to the main-flow velocity 4, the upstream 


disturbance can be expressed explicitly as 
S) = R. cos kihMiCA, A)r - (1x) + my) 
where (see Fig. 3) 


(4.1) 


l, = cos 6, -m, = snmé 


Eq. (4.1) clearly satisfies the governing equation for the 
entropy mode, Eq. (2.6). One notices here that a 
steady disturbance occurs as a special case at /; = 0. 

At a given orientation 6, the incoming disturbance 
drifts along the shock at a constant speed of 


 ] 


C, = [cos 6/cos (6 — €)|U; (4.2) 


Since the flow pattern of the upstream disturbance 
appears to be steady to an observer moving along the 
shock at this speed, the flow pattern downstream of the 
shock will also appear to be independent of time. 
With respect to the same observer, the downstream 
main flow, on the other hand, will appear to move at 


an apparent velocity of U’, defined by 


C= U0 +(=c) (4.3) 


Since the magnitude and the inclination a of LU, de- 
pends on the magnitude of C,, in particular, it is con- 
ceivable that U’, could be supersonic for some magni- 
tudes of C, but subsonic for other magnitudes of C,. 
One may thus introduce an effective Mach Number //, 


and a corresponding effective Mach angle yu, defined by 


M, = U,/A (4.4) 
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ue = arc sin (1/4,) (4.5) 


respectively. The dependence of .1/, on the shock 
strength x at a given shock angle 8 (e.g., 8 = 30°) and 
for a given inclination 6 of the upstream disturbance is 
plotted in Fig. 4. In view of these considerations, if 
one adopts a system of moving coordinates in conform- 
ity with the above-mentioned fact, the analytical work 
involved can be much simplified. This is particularly 
true if one selects a system of coordinate axes OX V in 
such a way that the X-axis and the Y-axis are taken 
along and normal to the velocity vector U’,, respectively 


(Fig. 3). (This scheme breaks down for upstream dis- 
turbance coming in a head-on direction—i.e., at 6 — 
« = 3/2. However, this special case can be singled 


out and studied separately.*: *) 

As has been mentioned previously, when an upstream 
disturbance of a given mode (as in the present case, the 
entropy mode) interacts with the shock, usually all 
three modes are produced at the shock. Since vor- 
ticity and entropy are particle-attached properties, with 
respect to the moving frame of reference OX Y, their 
trajectories, clearly, are along the velocity vector U’,. 
In fact, if one decomposes the downstream perturbed 
velocity into two components—# along OX and ¢ 
along OY, respectively —it can be shown that vortices 
simply arise from shear waves associated with the ve- 
locity component @, while the velocity component @ is 
purely irrotational. The sound mode or its equivalent, 
the scalar potential ¢, in the moving frame of reference 
of OX V is seen to be governed by the following differ- 


ential equation :*: * 


(M,? — l)odxx + dyy = 0 (4.6) 


Since Eq. (4.6) is hyperbolic or elliptic according to 
whether \/, > 1 or < 1, one reasonably would expect 
that the nature of the perturbed pressure field produced 
downstream will differ greatly in these two cases. Phys- 
ically, this implies that, for 17, > 1, the sound field 
generated downstream at a fixed field point is affected 
only by a localized distortion of the shock; while, in the 
case of IJ, < 1, it is affected by the whole shock con- 
figuration. In the former case, the pressure disturb- 
ance produced at the shock propagates downstream 
along a characteristic with constant amplitude (con- 
ventional sound wave). In the latter case, since part 
of the disturbance produced at the shock is fed back 
to the shock, its amplitude attenuates with the distance 
as it proceeds away from the shock. The essential 
difference between these two types of solution is illus 
trated in Fig. 6 for a normal shock. 

At this point, it might be proper to give a general 
outline for the steps to be taken in the analytical treat- 
ment; readers interested in the detailed procedures can 
refer to more extensive works done by the author.** 
\s has been just indicated, the nature of the governing 
equation for the scalar potential ¢ differs according to 
the range of the effective Mach Number ./,; conse- 
quently, our analytical treatment has to be divided 


into two main groups. However, instead of using the 
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PLANE SHOCK AND 


potential ¢, it has been found more convenient to use 
directly the pair of physical parameters @ and / (or its 


equivalent g = — cos 4, Pp) i.e., the two different as- 


q| _ | A cos p, S —m (MW cos 8B — C, 
é| | As: cos (a — B) | ~~* m3 (cos 8 — C 


at .Y sin (a — B) 


or alter the elimination of the shock angle yy, this 
gives a linear relationship between qg and 7% in connec- 
tion with the known function S i.e., the value of the 
upstream disturbance at the shock. 

Our general scheme is to use the above boundary 
conditions at the shock for the two flow parameters g 
and @ in combination with their governing equations to 
determine first the sound mode produced downstream. 
Then the corresponding shock angle can be evaluated 
algebraically from Eq. (4.7). Using the remaining 
boundary conditions at the shock, the value of the vor- 
ticity and entropy of the downstream perturbed flow 
field at the shock can then be determined. Their values 
in the downstream perturbed field now can be inferred 
directly from the governing equations. 

Specifically, when 1/7, > 1, g and @ are found to satisfy 
a pair of Wave equations; their solutions can be written 
down at once,’ thus 


Y cot u,) + FX + Yecoty,) | 


q = F(X = 
- Vcot u,) + F(X + Y cot u,)$ 


- — F(X 


(4.8 


t 


Since only one of the two functions F; or F: represents 
sound waves generated at the shock and propagating 
The 
selection of the proper one depends on the magnitude 
of C,—1.e., the drift speed of the upstream disturbance 
along the shock. It can be seen that, when C, is on the 
lower segment between the intercept of the shock and 
the sonic circle, F, is to be taken; when it is on the upper 


downstream, only one of them is to be retained. 


segment, / is to be taken (see Fig. 5). 

In the case of J, < 1 through a proper adjustment of 
the scale of the variables Y and Y corresponding to the 
usual Prandtl-Glauert transformation, the two param- 
eters g and @ are found to satisfy a pair of Cauchy 
Riemann equations. In conjunction with the boundary 
conditions at the shock, Eq. (4.7), the problem can be 
reduced to an integral equation of the Fredholm type 
involving one unknown function. With respect to 
the perturbed pressure field generated downstream, this 
unknown function may be interpreted as a specification 
for the distribution of dipole moments at the shock. 

Thus far, we have excluded purposely the special 
case of IJ, = 1. At J/, = 1, the governing equation 
for the potential @ reduces to 


dyy = OV (4.9) 


which is of the parabolic type. The two characteristics 
coalesce into one—1.e., XY = constant. 
wave the drift velocity is always normal to the wave 


Since in a sound 


front, this implies that 


OBLIQU 


A) 


+ 
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pects of the same physical entity, the sound mode. 

At the shock, after some algebraic manipulation, g 
and v satisfy the following system of equations: 


COS LL, , 
. sin (a — 8B) 
.1) cos (a — B) + mq sin (a — 8B) = vy ee 
(4.7 
VY cos (a — B) = 0 
v(or gy) = O and p = —F(x) 


Introducing these results into the boundary condition 
at the shock, Eq. (4.7), both the pressure perturbation 


and the shock angle ~y can be evaluated immediately. 


5) RESULTS 


For simplicity of illustration a normal shock is con- 
sidered. By taking the upstream Mach Number .\/, = 
1.45 (or the corresponding shock strength y = 2.286), 
a numerical example is worked out for an upstream 
entropy disturbance of the type specified by Eq. (4.1). 
As indicated in Figs. 7, 8, and 9, the amplitudes of the 
wave generated are functions of the inclination 6 of the 
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Fic. 11. Decay distance of outgoing sound waves with respect to 


orientation of incoming entropy mode. 


At \V/, < 1, the amplitude of 
the sound wave refers to its value immediately after 
the shock. 
Ay for the shock displacement is continuous at 1/7, = 1. 


incoming disturbance. 
The result indicates that the amplitude 


However, the slopes of the curve are not unique at that 
point; the right-hand slope is large but finite, and the 
left-hand slope is infinite in the neighborhood of J/, = 
1. A similar argument holds for amplitudes of the 
flow parameters themselves. A glance through the 
results obtained shows the following difference between 
the two types of solution : 

(a) When J/,> 1 the phase angles of the waves gen- 
erated downstream are in agreement with those of the 
incoming disturbance, while at 1/, < | there are defi- 
nite shifts in phase angles. 

(b) The sound waves generated downstream have 
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a permanent wave form when 17, > 1 but decay with 
distance when VW, < 1. Ata fixed value of the shock 


strength, the decay distance 
d = (1 — M?)/\Rkym [1 — .7,2]'7! 


clearly is a function of the inclination of the incoming 
disturbance. As indicated in Fig. 11, in general, the 
greater the inclination 6, the shorter the decay distance 
d. 

As a result of the interaction, the shock no longer 
remains a plane. Amplitudes of the corresponding 
shock displacement are shown in Fig. 10. 
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Inviscid Hypersonic Flow Past Blunt Bodies 


S. H. MASLEN* anp W. E. MOECKEL** 
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SUMMAR\ 


Two methods are shown for the calculation of the flow field 
etween a blunt body and the shock associated with it for the 
The 


symmetric 


Real gas effects are included 


that is, 


ise of hypersonic flow 


solutions consider only symmetric flows 
wdies at zero incidence 

One method consists in tracing successive stream tubes around 
the body and leads to iterations on the initially assumed position 
f the shock 
to the Karman-Pohlhausen procedure for boundary 


jistinction is made between round-nosed and flat-nosed bodies, 


The second is an integral method closely analogous 
lavers \ 


und both cases are discussed 
\ specific example corresponding to a re-entry missile situation 
is calculated; the two methods agree within a few per cent 


Comparison is also made with other known solutions in the 


stagnation region 


SYMBOLS 


parameters defined in Eq. (19) 
velocity profile function [Eq. (19) 
ratio of local to (constant) total enthalpy 
ratio of local to nose curvature 
density to stagnation point 


ratio of free-stream 


density, (p Po)x =0 
radius of flat-faced body 
V = Mach Number 
P = ratio of local to stagnation point pressure 
R = nose radius of curvature 

= ratio of local evlindrical radius to nose radius of 

curvature 

ratio of local to stagnation point temperature 
defined in Eq. (21) 
ratio of local velocities in a directions to free 
stream velocity 
coordinates (nondimensionalized — by 


the 


orthogonal 


dividing by nose radius of curvature), 


see Fig. 2 

shock-layer 
curvature 

a see Eqs. (19) and (20) 

= ratio of local to stagnation point density 


(47) 


bk = ratio of thickness to nose radius of 


f 


stream function [see Eq 


+ 


Subscripts 


value immediately downstream of the shock 


wall value 


partial derivative 
stagnation point 
= free-stream value 


2 


Jars denote dimensional variables 


INTRODUCTION 


- VERY HIGH VELOCITY FLOW past a blunt body, the 


shock lies close to the surface, tending to wrap itself 


Presented at the Aerodynamics—III Session, Twenty-Fifth 
Annual Meeting, IAS, New York, January 28-31, 1957. 
\eronautical Research Engineer 


** Assistant Chief, Propulsion Aerodynamics Division 
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In the calculation of skin friction 
to the body, the 
knowing the flow 


against the body. 


particularly, heat transfer 
boundary-layer theory 


field immediately outside the boundary layer. 


and, 
requires 


Several analyses of the external inviscid field have 


been performed under varying assumptions. Hayes! 
has discussed the problem in general terms and 
Li and Geiger*® have given a solution valid near 
the stagnation point. Serbin‘° has presented an 


approximate analysis for the subsonic region for a 
disc and for a sphere. Chester® has given a general 
solution of the problem for a perfect gas, as has Free 
man.’ The latter “boundary-layer” 


approach somewhat similar to one of the methods 


paper uses a 
presented herein. 

In the present report, two methods for determining 
the flow field about arbitrary body shapes will be 


described. Real gas effects are included. The solu 
tions consider only symmetric flows—-that is, sym- 
metric bodies at zero angle of attack. The first 


method is a graphical one, in which successive stream 
tubes are traced, and the second is an integral method 
analogous to the Karman-Pohlhausen procedure for 


boundary layers. 
THE STREAM TUBE METHOD 


The stream tube method was developed in reference 
S to compute, inviscid velocity and Mach Number 
profiles near the surface of blunted cones and wedges. 
The method is easily adaptable to the leading-edge or 
nose region at high Mach Numbers because the stream- 
lines are almost parallel to the surface except for a 
small region near the stagnation point. Furthermore, 
the surface pressures are known to be very nearly the 


Newtonian values (cf. Lees, reference 9). The pro- 
cedure can be outlined with the aid of Fig. |. We 
consider first the stream tube (owbc, Fig. 1) imme- 


diately adjacent to the body. This stream tube passes 
through a normal shock and so has a known entropy. 
With this and the known surface pressure distribution, 
the surface density, enthalpy and, finally, velocity, 
can be computed anywhere on the stream tube from 
data presented in thermodynamic charts based on 
real gas properties. Once the density and velocity 
are known, by continuity of mass in the stream tube 
owbc, the distance we is given by 


= (p./Bw)(u../ Mw) (ob? / 27 (1) 
where j, #, 7 are density, velocity, and cylindrical 
radius, respectively. 

Now, if the shock shape is estimated for the nose 
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region (say obd in Fig. 
continued for successive stream tubes. The entropy 
is known from the shock shape. The pressure can 
be taken as either the surface value or can be corrected 
for the centrifugal effect by 


P, = Kpu? (2) 


where A is the local curvature and 7 is the distance 
normal to the surface. With the known pressure and 
entropy, the remaining quantities can then be found 
and the process repeated until a point (e, Fig. 1) is 
reached where the mass flow across we is the same as 
that across obde. This gives the shock position at 
one radial angle 6. This process can be repeated at 
as many angular stations as are needed to determine 
the location and shape of the shock and the stream- 
lines. If the resulting shock angles associated with 
each stream tube differ from those assumed, the new 
values are used and the computation is repeated. An 
iteration on this basis will lead to a final determination 
of the shock and of the entire flow field. The principal 
source of error in this procedure is the fairing of the 
shock from the first station where the flow is nearly 
parallel to the body (point d in Fig. 1) to the axis 
(point 0). The location of the shock at the axis can 
be estimated quite accurately from existing theories 


> 


and experiments——e.g., reference 3. The fairing from 
o to d is best accomplished by making the shock as 
nearly parallel to the body as the fixed slopes at o and 
d permit. 

An example of a calculation made by this method 
will be presented later. For the present it suffices 
to say that only three iterations were required in the 


example in question. 


INTEGRAL METHOD 


Equation of Motion 


For the axisymmetric flow of an inviscid fluid, the 
equations of motion in an orthogonal curvilinear 
system can be written (this system is also used in 


reference 3) 


|), this procedure can be 
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(rpu)s + [(1 + Ky)(rpv) |; = 0 

uur + vu,(1 + Ky) + Kuv = —(P,/p) 

uv, + vv,(1 + Ky) — Ku? = [-( + Ky)P, | 5 “) 
h+ (1/2)(@? + 07) = hk, + (1/2)u,? 


where (see Fig. 2) *, ¥ are parallel and normal to the 
surface, respectively; K = K(x) is the local surface 
curvature; “ and @ are the respective ¥ and J velocity 
components; 7 is the cylindrical radius and P, 5, hh are 
pressure, density, and enthalpy. For plane flow, 
Eqs. (3) apply, but with 7 = 1. 

For an essentially round body, it is well known (e‘ 


Hayes!) that for high Mach Number flows, 
5/Ry = O(k) 


where 6 is the shock-layer thickness—1i.e., the distance 


between body and shock— R, is the radius of curvature 
of the body at the nose, and k is the (small) ratio of 
density upstream to that at a stagnation point behind 
the normal shock. The ratio so defined is related to 
the density ratio, k;, across the corresponding normal 


shock by 
k = ky [1 ai (ky 2) | + O(R,*) 
which result is found by considering that the gas 
behaves as a perfect fluid along the stagnation stream 
line between the shock and body. 
Similarly, for an essentially flat-faced body, Hayes 
has shown that 


6/L = O(Vk) 5) 


where L is, say, the radius of the face. 
Hence, in the manner of the boundary-layer theory, 
new dimensionless variables are introduced as follows 
u = ti/aU_, x = &/R 
v = 0/ki,, vy = ay/kRo 
p = p/p, = ee 6 
hk = P/P,, K = ad/kR, 
hk = Bijan * & = ER,/ a 


where a is | or Wk according to whether the body has 
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, rounded or flat face. The subscript ‘‘0” refers to 


stagnation point conditions. Under these transforma- 


tions, Eqs. (3) become exactly 
(rpu), + [(1 + Kky)(rpv) |, = 0 (7) 


plu, + vu,(1 + Aky) + Akuv| = [(—ke)/a?|P, (8) 





plkuv, + Rkvv,(1 + Aky) — Katu?| = 
—c(1 + Kky)P, (9) 
h+ou+ kv? =1+h (10) 
yhere 
o = Po/pi? = 1 — k/2 + O(R%, Mo) (1D) 


fq. (11) follows from the normal shock relations and 
the y-momentum equation |Eq. (9)| integrated, for 
isentropic flow, along the axis of symmetry. 

Further, if the Mach Number is very large and terms 
of order k® are neglected, Eq. (10) becomes 


h+ au? = | (12) 

Boundary Conditions 
On the body, one has simply v(x, 0) = 0. On the 
axis of symmetry u(o, y) = 0. From the symmetry 


of the flow it follows that P, p, h, and v are even in 
', while w is odd. At the shock, conditions are some- 
what more complicated. Because of the body-oriented 
coordinate system used, the respective normal and 


tangential velocities at the shock are not @, and “%,, 


but rather 7,, and #., (see Fig. 2). These are inter- 
related by 

U, = My COS € — 7, SIN € 

7, = d,, cos € + My Sin € 


From Fig. 2, 
€= tan~'[(db/d%)/(1 — K6)| = 
((d5/dx)/(1 — K6)| + 0(68/Ry*) (13) 
The oblique shock relations give 
Hga/u. = cos(¢ t+ e) 
da/u,. = [(—b.)/d;| sin(y + €) (14) 


Pipa? = [1 — (fig?/n,2)|[1 — (,/p.) | + OU /M?) 


lf Eqs. (6) are used, these last six relations yield 


u (cos g)/a| — (k/a?)6, sin g + O(R?/a*) (15) 
or 
“us = cos g — ké, sin ¢ + O(k?) (round nose) (15a) 
ue = [(r/2 — 9)/Vk| —6,+ 0(-Vk) (flatnose) (15b) 





P, = (1 — atu,2)(1 — k/p + k/2) + O(R?) (16) 


Equations for v,, analogous to Eqs. (15), can be written. 
However, this condition can be replaced by a more 
convenient integral condition. This is found by equat- 
ing the mass flow across AB (Fig. 2) to that across 


the shock in the region BC. The resulting expression 
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f = u(x, y)/us(x, 6) = a(x) + a(x) + o°c(x) 4 
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is, 1n dimensionless coordinates, 


(1/2)r,? = rpu dy (17) 
J0 

Eqs. (7), (8), (9), and (12) plus the boundary condi 
tions (15) and (17) and the condition at the surface 
and that of symmetry define the problem. As was 
pointed out earlier, two cases arise according to whether 
the body is flat-faced or not. The distinction depends 
on whether or not the difference between the shock and 
body angle is of the same order as the deviation of the 

shock from a normal one. 
It should be emphasized that the integral method, 
tube method, 


in contradistinction to the stream 


makes no assumption about the surface pressure 


distribution. Agreement between the calculated dis 


tribution and the known Newtonian value actually 


serves as a check on the method. 


Rounded Nose Solution 


-lnalysisIn this case a = |. It appears reasonable 
to try to solve this case by an integral procedure using 
a polynomial profile as in the Karman-Pohlhausen 
method for boundary layers. The solution to be pre 
sented in subsequent paragraphs is carried out to 
terms of order k, but, without any basic change in the 
scheme, could be arranged to retain terms of any 
desired order in k. 

First an integral of the x-momentum equation [Eq. 
(S)| is required. After some gymnastics, this can be 
written, correct to order k, as 


bd) 
| rpu-dy | — us|l + Kké| rpudy| + 
JV z J 0 Z 


%0 vA 


k[Kror — toK,] | ypu-dy = 2ku.u | rdy (1S) 
v7 


Next, let a polynomial profile be introduced. A cubic 
is used as being the lowest order equation which could 
exhibit the profile inflection that previous calculations 


by the stream tube method had indicated. Let 


a*d(x) 


(19) 


where a, b, c, d are functions of x (or, equivalently, 
u,), and where 


ey | 76 7 
¢ = | pr dy prdySy=t (da/pr) (20) 
« j J 


0 / 0 


t= pr dy S65 t (da/ pr (21) 
0 


If these relations are used in Eq. (1S) and u, is made 


the independent variable instead of x, then 


»! *1 
tu,” fda — u.[1 + KRd| | tu fda 


sd | *o 
(KK / Ter) ual? Ms” i | (da/p) | do 
. Ji 


2ku,.t (da/p (22) 
/ 
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Similarly, Eq. (17) becomes 
1 
tu, | fda = (f_"/2) + brat; (23) 
J 


where use has been made of Eq. (15a) and the results 
(from Fig. 2) that 


fur = +sing (24) 
/ » 
oe a Si ky Vv l — Teor" 
Indeed, from Eqs. (15a) and (24), 
Uy = V1 — Pur? — RKry,?bu, 
(25 
lor = V1 — uy? [1 — RKU,6,, | 


Also, it is convenient to rewrite Eq. (21) in the form 


“1 
§ = (t/rp) | (da/p) — (kus/2rw) X 


/7 
! 2 
| v3 } (do | (26) 
/J 0 


Finally, in the present notations, the y-momentum 
equation [Eq. (9) ] becomes, on using Eqs. (16) and (24) 
“1 


Piss | fPedo + O(R) (27) 


— u.> — (Ktuy?/ry) | 


Now, to determine the five functions a, ), c, d, and 6, 
that 
A third is that, from Eqs. 


five conditions are required. Two are Eqs. 
(22) and (23) be satisfied. 


(19) and (20), 


a+b+cHt+d=1 


(28) 


A fourth condition is on the surface velocity —i.e., 
on a. Because the streamline on the surface is that 


passing through the shock where it is normal, the 
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surface entropy is known. Hence, the enthalpy (anq 
from this, the velocity) is a known function of th 
surface pressure. To get from such a relation to om 
relating a and u,, we proceed as follows. From Eqs 
27), (22), and (23), there follows successively, 


.77) *1 
P, =1—a; — K/s,) E : } Pao | du, + O(k 
0 70 ts 


mu | 
= 1 — u,° — (K/r, | uy E fda x 
J0 J 0 ju 
du, + O(k 
= | — u,? — (K/r,) | Ul yV ingly + OCR) (29 
J 0 


If terms of order k were retained at this point, the 
they should have been retained earlier [in Eq. (27 

Such effects can readily be found by iteration but will 
not be needed at 


present. This point is discussed 


further in Appendix A. At any rate, Eq. (29) and the 
known wall temperature-enthalpy relation determines 
the surface velocity (a@ versus u,). 

The final boundary condition is one which determines 
the wall velocity gradient in terms of the shock curva- 


ture. In Appendix B, it is shown that 


uy(x,O) = putwly [((1 — 2k)(Ro/R;)? — 
k(Ku/rpT),-0| + O(R?) (30 

where the R,/R, is the ratio of the radii of curvature 

of the body and shock at the stagnation line and is 

given in Eq. (B-4). Using Eqs. (19) and (20), this 

leads to 

_ Dp 


Gg: / th.) Ves 


1+ 6+ 


-k(Ku prT|,-04 (31 


68, /t lo = Girehil 





If Eqs. (22) and (25) are combined with Eq. (31), 


one can write, for later convenience, 


*) ‘ 
3b { Pde = Fy=(3T, us?) }1 — 2k[1 +6 4+ (6, Us) uso — R[Ku prT|,-0t X 


| Mata us(l + KRb) + uR(6U re) ve + 
J 0 


>! 
2b | fdo — F, => CP ofan” 


If the integration on the left-hand side of these expres- 

sions is carried out, the results can be combined with 

Eq. (28) to yield 

b4 — (l5a + 3)b* + (72a? 
[(66a + 12)F% 


— 12a + 18 — 6F,)b? — 

+ TOP; |b + 72 Fe? = (34) 
Once 6 is known, the remaining quantities can be found 
conveniently by solving Eq. (31) for ¢/r,, and then, 
from Eqs. (23) and (28), 


~ $orw/tUs Ou, 
¢ = 3° 1 + 2k — 1|1—b— 8a> 
' t Vu Vu f 


—a—b-—-ce (35) 
The necessary equations for the solution of the flow 
about a rounded blunt body are now set up. 


U7), 1 + 2k(u6 ry) — 2R{1 


“1 * 
QRriptsd + RK ty) ug Pus? } E } (do | do, du, (32 
( 7 


J0 


+ 6 + (6,5 Us) |p-0 — R[ Ku ‘prT |,~05 (338 


The Solution Procedure--A practical method for 
solving the equations follows. 

(1) Compute a from Eq. (29) [using Eq. (25)] and 
the known wall pressure-enthalpy relation. Also 
find 7), and py. 

(2) Now consider Eqs. (34), (31), and (35). Ii F 


and F, and an approximate value for 6 are known, thes¢ 


expressions yield b, ¢, c, and d. The problem is to 


get a good estimate for Fi, /2, and 6. Everything 1s 


known in F, and Fy [Eqs. (32) and (33)] except 6 


6(u;) and | (da/p). Both of these are small 


0 
terms (proportional to k). For 6, some experience 


shows that a good initial estimate is 
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— H,”) (36) 


where 6 is the value of 6 at the stagnation line and can 
be taken from reference | or 3 or the stagnation region 
analysis which follows this section. 

The remaining term to be estimated is 


*1 *o 
I, = WKB / te) ae f°," ii (da/p) | da 
J0L /J0 
This is a very small term so a rather rough estimate 
will suffice. Thus, 





as BO) Folas te 4 (4/ Ha) (da/p)| X 
e ( 


) 


>] 


~ RK / ta) ue 6,4 | (t/U,) | af-da 
J0 


bd | 

The quantity ¢/u af-do is fairly constant, in- 
JW 

with uy. <A 

approximation is to use its (constant) stagnation region 


FF, and Fs so 


c, and d can be found from Eqs. (54), 


creasing slowly sufliciently accurate 


value. This determines /; and hence 
that b, t/rip, 
31), and (35). 

(3) Next, to find the shock-layer thickness 6 [Eq. 
This 


and the pressure. 


26)|, the local density variation is required. 





follows from the enthalpy [Eq. (12) 
There is no difficulty with the enthalpy. However, 
the pressure is another story indeed. The wall pressure 
has been used in the form given by Eq. (29) which 
neglects terms of order k. The present determination 
of f includes effects of order k. These tend to increase 


the magnitude of the term [see Eq. (27) | 


*| 
— (Ktu,?/r,) fda (37) 


However, to order k, there should be other terms in 
which will tend to 


27)—see Eqs. (9) and (16) 


Eq. 
counteract this effect. Thus, the shape of the variation 
of pressure across the shock layer is still given approxi- 
mately by Eq. (37). Hence we write, as the pressure 
variation in the shock layer, 


P= P (A /2r,-) | Ue? — fe a 


fda | rac) (3S) 


Where it will now be assumed that this equation is 


correct to order k. This assumption must be checked 


at the completion of the calculation. Using this 


equation and Eq. (12), P and / can be found and, from 
them, p. 

(4) Next 6 is given by Eq. 
be compared with Eq. (36) and the solution iterated 


(26). This result is to 
if necessary. 
(5) Finally y is related to o by Eqs. (20) which can 


be written, analogously with Eq. (26), as 
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*o 
y =] (t/re) (da/p) | — (kuy./2ry) X 


70 


LZ / Tes} (da p) (39 
o 


The approximation in the pressure [Eq. (3S)] can 
be checked to complete the solution. This is done by 
integrating Eq. (9) across the shock layer, retaining 
terms of order k and by retaining terms of order k 
in the shock pressure [Eq. (16)|. If the resulting 
surface values agree with those found from Eq. (29 
the matter can be considered closed and the problem 
solved. For the example to be discussed such is the 
case. Actually, the surface pressure distribution is 
known to be Newtonian. Hence, if Eq. (29) indicates 
such a variation, there is really no great need to pursue 
the check just described. Of course, if one carried out 
the present analysis neglecting terms of order k every 
where, then Eqs. (24), (25), and (29) yield the New- 
tonian pressure only if the integral term in Eq. (29 
is omitted. However, for calculation, an approximate 
relation [Eq. (38)| for the pressure has been used as 
A check on this 


admittedly crude assumption is provided by comparing 


though it were correct to order k. 
Eq. (29) with the Newtonian values. The agreement 
which occurs arises because the pressure reduction 
due to consideration of centrifugal effects is very nearly 
balanced by the increased pressure due to the shock 
being stronger than if it were parallel to the body 
surface. 

Stagnation Region—Betore proceeding with a specific 
example, it is of interest to examine the present solution 
in the neighborhood of the stagnation region where 
several other known solutions (Hayes! and Li and 
Geiger*) exist. In the present analysis, the neglect 
of various terms depends on distances normal to the 
surface being substantially less than those parallel, 
and hence that 0 OV > O Of. 


in the stagnation region. However, 


This is not the case 
the stagnation 
region is not really singular (as contrasted with the 
leading edge in a boundary-layer flow) and the solu 
tions can be written in a Taylor series (in x) about 
r = 0. 
basis yields the information that the terms heretofore 


A formal expansion of the equations on this 


neglected are still the negligible ones near the stagna 
tion region. A check on this is provided by comparing 
results of the present analysis with those of references 
| and 2. 

In the stagnation region, the boundary condition 
determining the wall velocity is particularly simple. 
If one applies the x-momentum equation |[Eq. (S) | 
at the wall, there follows successively, using Eqs. 


(19) and (29), 


Py(Ua) (Ua), = 


or, as x — O, 
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Fic. 38. Shock-layer thickness on stagnation line 


a> V(S 3k (40) 


Next, in this region, /, and F, |Eqs. (32) and (33) | 


can be simplified. First, using Eqs. (25) and (B-3), 





Uy > ry[l — (Rb,/x) (41) 
then /; and Fs become 
*! 
Fi =1—k 2 +a — 36+ & | ve | 
J0 r=0 (42) 
Fo = 1—k(2+ a) 


To evaluate these quantities, an initial guess for 6 and 


“1 

for | afdo is needed. The former (6) can be 
J 

taken from references | or 3, and the integral is about 

1/4. and then 3, t 7, 

c, and d [Eqs. (34), (31), and (35)| can be found. 


Using these values, 7; and F 


Then 6 and y are found from Eqs. (26) and (39), which 
for this case become simply 


= (te) (R/2)(t/tr»)? 
(45) 
(k/2)(to/r»)? 


v= (to ry) 


The variation of shock-layer thickness is shown in 
Fig. 3 and compared with Li and Geiger’s result.’ The 
Note that to this order 
the results are independent of the nose shape. 


differences are rather small. 


Fig. 4 shows the variation of the profile function 
f [Eq. (19)] for k = 0.073. 
with that of Li and Geiger. 


The result is compared 
The 
are mainly due to the curvature effect on the pressure 


small differences 


Li and Geiger have a = +/ 2k, because they neglect 





the integral term in Eq. (40). | 


Flat-Faced Bodies 


In describing the general problem of the flow over 


blunt bodies a distinction was made earlier between 


rounded noses and flat-faced ones. A rough analysis 
for the flat-faced case follows. We consider here only 
the first-order solution where k > 0. Then Eqs. (9), 


(10), and (16) yield (a = Vk) 


P=h= p= 1+ 0(k) (44) 

Then Eqs. (7) and (8) become 
(72), + (rv), = 0 . 
(45) 


uu, + vu, = 0 


Also, because the curvature is limited to be at most 
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of order /k [Eqs. (6) ], it follows that 


tf = 


x — kg(x) 


K = (V 2z¢,), 


and 2/L = Vk | V 29, dx 
e ( 


) 
where 2 is the streamwise ordinate of the body surface 
From the first of these equations it follows that ¢() 


must vary at least as x* for small x. Then, to th 


present order, x can replace r in Eqs. (45). If a strean 


function, w, is defined such that 


xu = yp, : 
xv = —y, ii 


the general solution of Eqs. (45) is 
Y = vlxyy + fi(x) | 


fi(x) is arbitrary but is zero if it is required that v u—>() 


as y— 0. Hence, 
y = (ry) 
uw = P'(xy) (48 
v = (—y x)v' (xy) 
Now the shock conditions are to be invoked. From 
Eqs. (15b) and (17), using Eq. (46), 
y'(xb) = V 2¢, — 6, 
(49 
o(xb) = x7/2 
Eliminating y, these can be combined to yield 
x + (6 + x6,)(6, — V 2g,) = 0 (90) 
For g, Q, this is the same result as that of Hayes 
Then 
F “ meee fa 3/4 
x = —6,[(1 + 26,”)/3 ; 
a : — P { 1 
6 = (1 + 6,7) (C1 + 26,7)/3] 


where the boundary condition used is that 6,, be per 
Then 


If g, does not vanish this becomes 


mitted to become infinite at the corner (x = 1). 


As => i) = —1. 


6, = ¥V 29, — | 
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Fic. 4. Stagnation line profile function for k = 0.075 
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‘ If Eqs. (49) are used, this can be written 
@ 
¢ P,—- 1 2 7 
| ' A : = —(V2¢, — 8,)? 4 x6 X 
| — a el k x" J0 
ee Be RL 
| V (Yo 
‘ -> 
4} | (V 25 6,)dx x 
| Cc \ 
BODY SHAPES "a 
(NV Dg 6,)ax ? 18) 
. Jt 
riace . 
t g 3 ; A For the case where Eqs. (53) and (54) hold, this yields 
x CONCAVE ; 
O the al~ | » 2 9 21/2 7 
B (Po — 1)/8 — (4x7/3) + 10 1) 6] * Gil 
rean . 2 
w FLAT (1 1) G[Ge(x) + [1 l G3 (A oY 
= 
Y 
n , 
o where 
7 a. 4 
‘ 
w | G(x) fIn{[(1 + y)/ 2] (3/°2)x* sin-' x 4 
< C 
a CONVEX 2x fo(X) (1/30)|—30y° + 27 + 30y' 
a ~80 2 4 6 8 10 5 (ee a 
; IOV O0Oy l6ov + 148 
x= K/L - 
5. Surface pressure distributions for several almost flat G(x) 14 In|(1 Tr Y) ? 21(1 y (1 - y 
ay faced bodies ' - 4 : 
v8fo2(x) + (1/15) [— ay Sy? + 120y! +4 
lOv 2394 LoOy 52 
c . , 1 2 rye . ° - 
\lso 6,, diverges as (1 x . Then it appears 
; : oa : - : ; ' 26 InI( | /) 9] Oy : ~ y 
(48) | that for a suitable variation of g(x), 6 will have a series Gs 6 In| C1 yije, vr se" sims 1 5) X 
eae — [10v8 — 9y® — 75y! — 30v* + 210y? Bay + 29 
lution in powers of V | x? where the coefficient 105 : ih 10: 15099 ‘ 
rom | of the term proportional to V 1 x? vanishes. If and where 
\ 2g x, the solution of Eqs. (50) and (52) is 6 = 1. 
- as . . y Vv i x*,  fo(x) y(3 + 2y)/(1 + y) 
For other variations the solution is not apparent. ia 
19 However, in Appendix C an approximate solution 1s ae © theer become 
given for the case 
Gi(x) —> —7/8 + (ix?/24 
+) = 4) . le E ‘ , : 
V 22: I - Go(x) > 77 16 — 83x? 96 
5() Then G(x) —> -—-4 S4+ x° 2 
) 
‘ =2-—.1+ [1 —.1) 2]d ¥?) - In Fig. 5, there are plotted suriace pressure distribu 
, : - tions for bodies of several different shapes. The point 
V (1 - 1) 6[( - \ l+ Vl-—-wx (54) : ee aA : 
of consequence is simply that dishing the face decreases 
' For . 1 | this is exact and for .! = 0 it agrees with the magnitude of the pressure gradient. In this case, 
Eq. (91) to within less than one per cent. Then from the solution has not been carried to high enough order 
r- | Eq. (49), ¥ is given by (terms ol order Yk are neglected) to give more than 
- qualitative results. 


pry = V2¥yj)2 — A + [C1 — A)/2](1 — 2) - 
| Vii — A) 6[U — 2y)?" a + V1 — 2y)if (55) 
| Using these results, the surface pressure distribution 
to order k can be found. Eq. (16) yields, noting Eq. 
4), Pp, = 1 — bu,? — k/2, 


7), (48), and (46) there follows successively 


so that, from Eqs. (9), 


| 


| 
| fy = 


| Pp, = 1 —k (V 2g, — 6,)? + (2k x?) | 


Ku? | 





— pk|uv, + vv, — 


bd) 


yy "dy = 
kk | y""dy — k/2 
/0 


— k(V 2g, — 6,)? + (2k/x?) xyy'dy — 
R('V 2¢,),/x| } y'dy — k/2 
0 


EXAMPLE 


The flow field has been calculated by both the stream 
tube and integral methods for an hemispherical nose 
moving at 25,000 ft. per sec. at 100,000 ft. altitude. 
The static pressure, temperature, and density ahead 


of the shock are 


P 0.0106 atm. 
T 392°R. 
B.. = 0.001066 Ib. /ft.* 


The calculations were made for air in equilibrium 
using thermodynamic charts based on the computed 
properties of air presented in reference 10. For the 
stream tube method, the exact oblique-shock relations 
of reference 11 were employed. Using these references, 
the density ratio and the stagnation pressure and 


temperature are 
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Fic. 6. Streamlines and shock calculated by stream tube 


method. 
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Fic. 7. Shock-layer velocity profiles 
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The results are shown in Figs. 6 through 9. Fig. 6 
shows the streamlines and shock location obtained by 
the stream tube method. The dashed portion of the 
shock is the faired portion. The streamlines are 
seen to be almost parallel to the body for values of @ 
greater than about 25° which was the initial station 
for these computations. 

Fig. 7 shows velocity profiles for several stations. 
[The main property of interest is the approximate 
linearity of the profiles especially near the wall. In 
Fig. 8 the profiles are again shown for several stations 
and the results compared for the two methods of 
calculation which have been described in the present 
report. The differences in velocity are quite small 
throughout. 

The variation with x (or @) of several characteristics 
of the shock layer is shown in Fig. 9. For comparison, 
the shock-layer thickness as calculated by the two 
methods is shown. The two methods give virtually 
the same result up to about 45°, where the difference 
is about 2 per cent. The surface pressure depicted is 
that calculated by the integral procedure and agrees 
with the Newtonian result (? ~ cos? @) except at the 
larger angles (6 per cent deviation at 6 = 60°). This in- 
dicates that the pressure reduction on the surface is just 
about cancelled by the increased pressure due to the 
shock being stronger than if it were parallel to the sur- 
face. This result is in agreement with results quoted by 
Lees.’ Finally, the surface velocity and Mach Number 
ie, for this case at least, virtually linear over the 
range of calculations. The results shown for these 
quantities and the sonic line were all obtained by the 
integral method. However, the two methods gave 
results which in all respects agreed to within a few 
per cent. The sonic point on the body is at about 
2° from the axis, while that on the shock is at 18S 
rhis illustrates the very high velocity gradient across 
the shock layer. For example, at 18°, where the 
Mach Number behind the shock is 1.0, the surface 
value of the Mach Number is 0.4. This change occurs 
over a distance of | in. for the present example if the 
nose is considered to be 3 ft. in diameter. Although 
this is extremely large shear when thought of in terms 
of inviscid flows, this is not the case when a boundary 
layer is considered. In the present example, the 
inviscid velocity change is about equal to that across 
the boundary layer. However, the ratio of boundary- 
layer displacement thickness to shock-layer thickness 
is for the present case (Ry = 1.5 ft.) of order 10~* so 
that this shear outside the boundary layer probably 
has little effect on skin friction and heat transfer. 

Although the shock will remain quite close to the 
body for larger angles than are indicated in Fig. 7, 
the methods do not appear to yield reliable results. 
In the case of the integral method, e (Fig. 2) becomes 
large and the approximations indicated by Eqs. (13), 


(15), and (16) become poor. For the stream tube 
method, the three iterations used were insufficient to 
obtain final shock shape for @ greater than about 60 

However, since this value of @ is already well into the 
supersonic range, the solution could be continued by 


the method of characteristics. 
DISCUSSION 


Two methods have been presented for the calculation 
of flow fields in the region between a blunt body and 
its accompanying shock for very high-speed flow. The 
two analyses yield results which agree within the limits 
of accuracy expected for either procedure. In both 
cases the calculation effort required is moderate even 
when real gas effects are admitted. Actually real 
gas effects do not increase the difficulty of solution 
over that for a perfect gas. However, the present 
calculations were made on the basis of an equilibrium 
state. If recombination rate equations were to be 
considered the process would be considerably com- 
plicated, although such a generalization seems feasible 
since the processes, in both cases, essentially involve 
only forward integration from the stagnation region. 

It is this question of forward integration that intro- 
duces a fundamental limit on the integral analysis. 
The flow equations have been treated as parabolic 
rather than, in the subsonic region, elliptic. This 
is valid so long as terms of order k”’~ are negligible 
(see Appendix A) compared with unity. If such is 
not the case, another boundary condition depending 
on downstream conditions is needed and the integral 
method fails. In this connection, if the present cal- 
culation is continued far around the body, the results 
become less and less reliable. This is because of both 
the above-mentioned dependence on a downstream 
condition and the failure of some of the approximations 
as pointed out in the example. 

In the stream tube method this limitation does not 
arise directly. The downstream condition is replaced 
by specifying the surface pressure (the modified 
Newtonian value). The calculation can readily be 
made accurate to the order to which this distribution 
is valid. As is pointed out by Lees’ the Newtonian 
distribution is quite accurate unless the local surface 


inclination is small. 
APPENDIX A—-THE y-MOMENTUM EQUATION 


If terms of order k& are retained in the y-momentum 
equation [Eq. (9)], that equation becomes (a 1) 


P, — pKu? = —kp[uv, + vv, + K°u*y (Ku 


One can eliminate v from this equation by the use 
of the continuity of mass expression [Eq. (7)]. The 
result is not particularly inspiring to behold, but has 
the important feature that it contains second deriva 
tives with respect to x. This unfortunate occurrence 
causes a basic difficulty because it then means that 
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the final integral equations are really of second order 
in x (or u,) and, hence, that forward integration will 
not work. This follows because there is only one 
available condition at the axis of symmetry. 

What all this amounts to, of course, is that the 
subsonic flow is really elliptic, not parabolic as has 
really been considered here, and the sonic line cannot 
really be determined independently of the downstream 
(supersonic) region. However, such terms as raise 
the order of the equations are of order k in Eq. (A-1) 
and will introduce terms of order k*’” in a [Eq. (19) ]. 
[This is most apparent from consideration of Eq. 
(40), although it follows generally from Eq. (S).| 
Hence, in the present application the error is of order 
k®* and the retention of terms of order k elsewhere 
(in the x-momentum equation) is justified. 

A parenthetical remark on the implications of this 
discussion with reference to the stream tube method is 
perhaps in order. Although, again, the flow really 
depends on the downstream region, such dependence 
is replaced in the stream tube method, by the specifica- 
tion of the surface pressure (the modified Newtonian 
value). Then the accuracy of the solution is limited 
only by the order of exactness of this distribution. 
This accuracy appears to be very good so long as the 
local surface inclination is not too small.’ 

Now let us return to the y-momentum equation. 
With reference to the discard of terms of order greater 
than k, it might be pointed out that the same (and 
more) terms are neglected in Li and Geiger’s analysis 


of reference 3. If, in their notation, one writes 


u = xf'(y) + x8g’(y) +... 
then their Eq. (85) becomes 
f''( + Ky)? + 8[g’ + (f'2/f)] + 2kf? + ... =0 


where the g’ term is the largest of the terms arising out 
of (uv, + vv,) in the y-momentum equation and 2K/” 
is the largest curvature term. If the derivatives are 
made of unit order by changing from y to Ay/k, then 


this becomes 
f'’'(1 + R(Ky/k)|? + (S8k?g’/K?) + 2kf" +... = 0 


so that the g-terms are small (of order k*) but the 
curvature term (2k/”) should be 
essentially the retention of this term in the present 
analysis that leads to such differences as appear in the 


retained. It is 


stagnation region comparison shown in Figs. 3 and 4. 


THE NORMAL GRADIENT OF THE 
SURFACE VELOCITY 


APPENDIX B 


A complete expression for (u,),~9 will be found (for 
a = 1), although such exactness is not used in the 
present analysis. The change in entropy, in terms of 
pressure and enthalpy is, for a real gas 


dS = (1/T)|[dh — (dP/p)| 


Because the entropy is constant along streamlines, this 


leads to (see Fig. B-1), 
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Na — he — (2he1/p.-)(P. — P-) = (Te/Tu) X 
[he — hy — (2hes/pu)(Po — P. 


But, from Eqs. (10), (14), and (9), 


Na — he = —(AB)*41 — (B?/p.2)[1 + (p95 p-)]! 
P, — P. = —((AB)?/alil — (Rp) [1 + (pss 2p 
hy — hy, = —2u,Au 


Py — Py = (Kpw/t1)Uy*Ay 
Also, by continuity of mass through the stream tube 
cabw, 

(AB)? = 2p ytty( Ro R,)°Ay (B-] 


where R, R, is the ratio of the radii of curvature o} 


body and shock at the nose. If these last six equations 


are combined, there follows 
(Uy)o = (Corl) Teli {1 — (Rp) |2(Ro/ Ro)? 
k(Ku prT),7 4 (B-2 


The ratio Ry R, is needed next. At the nose {see 


Fig. (B-1)] 
R/R, = 1 + 7,2)? /7,, = (—W 1 — Py?/tyy 


where s is the distance in the streamwise direction and 
x’ is the distance along the shock. Near the nose 
dx’ = (1 + kK6)dx 

Also, by virtue of the nondimensionalization used, 
To = xX — (x2/3!) + (25/5!) (d*7,,/dx*)-~0 + (B-3 
Finally, from Eq. (24), 

r= ty + kROV 1 — Pus? 
The last five equations yield 


R./Ry = {C1 + kd)2/[1 + k6 — (B6,/x)|},-0 (Be 


Equations (B-2) and (B-4+) give the complete expres- 
sion for u, at the wall in terms of the shock-layer 
thickness (6), its curvature (6,/x), the shock strength 
(k) and the geometry and local wall conditions. For 
the present purpose, the result will be used only to 
order k. 

From Eqs. (19) and (20) at the wall, 


(ty) = (4:0/b) Ore 
which, when combined with Eqs. (B-2) and (B-4) yields 


b = (t7,,/us)} 1 — 2k[1 +6 + (6,,/ts) uso — 
k{Ku/ prT |, of + o(k?) (B-5 
For the two-dimensional case, Eq. (B-1) is replaced by 
AB = pul wtty(Ro R,) Ay 
from which it follows that Eq. (B-2) is replaced by 
(ty) = —R(Ku), 


so that Eq. (B-5) becomes simply 


b = (kt/u,)(Ku/p),, 
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\ppENDIX C-—-APPROXIMATE SOLUTION OF Eg. (50) 


If Eq. (53) holds, Eqs. (50) and (52) become 


x + (6 + v6,)(6, — Ax) = 0 (C-1) 
$(x = 1)=A-—1 i 
“ (C-Z) 
ie = 1) =2=— 4 
\lso, 6,, diverges as (1 — x2)~'’*. Then it appears 


that 6 has a power series of the form 


6= > a,(V1 — x)" (C-3) 


=) 


Variation with x* rather than x is used because the 
glution will then have, explicitly, the desired sym- 
metric form with respect to the stagnation line, x = 0. 
If Eq. (C-3) is put into (C-1) and Eqs. (C-2) are noted, 
the first few coefficients in Eq. (C-3) are found to be 


mH=2+A, a =O, a = (1 — A)/2 


a, = —V(1 — A)/6, a = (17 — 9A)/48 
— =(112 = 17 + 8112/5) 
_ 14 V6(1 — -1) 
527 — 2794 «112 — 117A + 8142/5 
COS 7776(1 — A) 


Numerical values of these a, are given below for several 


values of A: 








A=0 A= -1'2 A=-+1/2 

a 2 2.50 1.50 
ay 0 0 0 

a: +0.500 +0.750 +0.250 
ds —0.410 —(0.500 —0.289 
ay +0.354 +0.448 + 0.260 
(ls —0.319 —0.405 —0.231 
(lg +0.291 +0.371 +0.209 





Because the coefficients are alternating and varying 
slowly in a consistent manner, it is suggested that one 


would write 


b= a) + ao(1 — x2) — a3(1 — x2)? "11 —-V1l-x+ 
1 — x? — (1 — x)? * + 


so that 


g=2-— 44+ (1 — A) 2] — x?) - 


Vl — A)/6f(1 — x?)?/7/[1 + V1 — x2]} (C-4) 


This gives the exact answer when A = 1(6 1) and 
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Shock 





Stream- 





line- = 
| a [ &- 
—~-— J 
a ” 
1<—— x, ———1m 
Fic. B-1 
can be compared with Eq. (51) when .! = 0. Thus, 
Eq. (51) Eq. (C-4) 
6 (0 2.28 2.29 
6, (O) 0 0 
5,, (O) —0.478 —().488 
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Flutter of a Helicopter Rotor Rotating 
in Its Own Wake’ 


R. TIMMAN* anp A. I. vAN bE VOOREN** 


National Aeronautical Research Institute, Amsterdam 


SUMMARY 


A theory is presented which enables the calculation of the 


aerodynamic forces acting on a rotor which is rotating in still air 
It follows that the usual formulas for an oscillating airfoil in two 
dimensional flow may be retained, provided a new circulation 
Flutter that 
regions of flutter become possible. The flutter frequency is an 


function is introduced. calculations show new 


approximate multiple of the rotor speed, which means that an 
accumulation of vortices in the wake takes place. Experiments 
with a small two-bladed rotor model confirm the new theory 
Jecause of the low Reynolds Number, this model shows some 
additional regions of instability. These can be explained by the 
known reduction in aerodynamic damping which is due to the 
chordwise motion of either the point of boundary-layer transition 
or the point of laminar separation. These additional regions of 
instability vanish if turbulence wires are fixed to the model. 


SYMBOLS 


g = structural damping 
2Al = 


distance between the mid-chord points of successive 


blades 


2irk/N = phase difference between two successive blades 

kp = stiffness against flapping 

ky = stiffness against pitching 

! = semichord 

} = distance of chordwise section to rotor axis 

t = time 

w = downwash 

x = chordwise coordinate 

8s = product of inertia of the blade with regard to 
pitching and flapping axes 

I = moment of inertia of the blade with regard to flap- 
ping axis 

Pi = moment of inertia of the blade with regard to 
pitching axis 

Liz = aerodynamic forces 

N = number of blades 

| = circulation function 

R = radius of the rotor 

a = flapping angle 

y*(x, t) = bound vorticity 

e*(x, 1) = free vorticity 

él = distance of pitching axis aft of quarterchord axis 

¢ = pitching angle 

A = (2/N)[(k + (»/Q)] 
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v = frequency of oscillation 

VR = natural frequency in flapping 
V7 = natural frequency in pitching 
p = air density 

w = reduced frequency, v//Qr 

Q = angular rotor speed 


(1) INTRODUCTION 


Ir AN EXPERIMENTAL INVESTIGATION of the dynami 

the Kolibr; 
helicopter, it was found that unstable oscillations could 
Th 


behavior of a model of the rotor of 


develop if the model was rotating in still air. 


Kolibri (Fig. 1) is a two-bladed ram-jet helicopter 


built by the Netherlands Helicopter Industry, whicl 
has the special feature of a relatively weak elasti 
connection in pitching between each blade and th 
on the model 


rotor hub. The oscillations observed 


had nothing to do with ground resonance, since the 


rotor axis was completely fixed, while, moreover 
stroboscopic illumination revealed that the blade 
oscillation was perpendicular to the rotor plane 
Clearly the oscillations were of the flutter type. The 


flutter frequency appeared to be always near to a 
multiple of the rotor angular velocity. If this multiple 
was even (2 or 4), the oscillations of the two blades 
were in phase, while, for an odd multiple (3), the two 
blades were in antiphase. This means that, at any 
point of the rotor circumference, each blade performed 
the same motion and generated the same vortex pat- 
tern. Hence, it was concluded that accumulation o! 
vortices in the wake was an essential condition for the 
occurrence of flutter. This conclusion was confirmed 
by the fact that a small air speed perpendicular or 
parallel to the rotor plane was sufficient to prevent the 
instability. Consequently, this type of flutter will 
only cause practical difficulties at the ground or in 
conditions of flight near autorotation. 

Calculations emploving the usual two-dimensional 
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aerodynamic coefficients did not show the observed 
instability, which is not surprising since these co- 
eficients will be modified considerably by the process 
of vortex accumulation. 

The first author (see reference 1) has previously 
established a two-dimensional theory for the aero- 
dynamic coefficients of a rotating rotor by considering 
the flow on a cylindrical surface of which the axis 
coincides with the rotor axis. Unrolling this cylindrical 
surface on a plane, a cascade of airfoils is obtained, and 
the corresponding two-dimensional flow was calculated 
with the aid of a conformal mapping method for the 
case that all airfoils (rotor blades) move in phase. 
Recently a more simplified theory has been given,’ 
which is briefly explained in this paper. In the two- 
dimensional picture of the flow the integral equation 
jor the distribution of the bound vortices on the blades 
and the free vortices in the wakes between the blades is 
sven. Neglecting mutual interference of the biades, 
this integral equation is simplified to an equation of the 
same type as the equation for a single oscillating airfoil; 
the only difference being that the blade now has a wake 
of free vortices upstream and downstream. The Kutta 
condition then yields a relation between the intensities 
of successive free vortex sheets. Since, in principle, 
the equation has the same form as in the classical case, 
The mathematical 
without 


its solution follows the same line. 
formulas are given in Sections (2) and (3) 
details of derivation (these details are given in reference 
2). The result of this theory is that the usual formulas 





for the coefficients for a single airfoil may be used, but 
the influence of the two sheets of free vortices necessi- 
tates the substitution of a modified circulation function 
F. 

In the complex P-plane this modified function travels 
along the circle with radius 0.5 about the point P = 0.5. 
This means an important modification of the circulation 
function since, for a single wing, the real part of P is 
always larger than 0.5 and the imaginary part always 
negative. 

The three-dimensional configuration of the wake has 
been taken into account by Greidanus, Timman, and 
Zaat® who obtained the result that, for rotors with the 
usual small solidity factor, the two-dimensional approxi- 
mation is satisfactory. 

Calculations‘ confirmed that the modification in the 
Pfunction has an important effect on the flutter 
behavior. The rotor investigated had two stiff blades 
which could perform a flapping and a pitching motion. 
The inertia axis coincided with the pitching axis at 20 
per cent of the chord aft of the leading edge. If the 
natural frequency in flapping was larger than that in 
pitching, the system was found to flutter in certain 
angular speed ranges, according to both calculations and 
experiments. However, for flapping frequency smaller 
than pitching frequency, the system should theoreti- 
cally be stable, but the experiments showed that flutter 
still was possible in a mode which chiefly consisted of 
rotation about the pitching axis, thus suggesting flutter 
of one degree of freedom. 





695 





WAKE BLADE 


-2h-1 -1 +f | 
_Dhel 2Qhel 


Fic. 2. Unrolled cylindrical surface with blades and wakes 


The theoretical aerodynamic forces appeared to be 
such that, although flutter with one degree of freedom 
is impossible, its damping is much reduced compared 
with the case of a single wing. The theoretical forces, 
for boundary-layer 
Number 


however, do not account 
effects. At the relatively 
(~32,000) for which the one degree of freedom flutter 
In reference 


any 
low Reynolds 
occurs, such effects may be of importance. 
5 it was shown that, in the region of critical Reynolds 
Numbers, important differences between the experi- 
mental and the theoretical values of the moment about 
axes near the quarter-chord axis are possible. In 
particular, the experimental coefficients showed lower 
values for the damping. As will be explained in more 
detail in Section (4.4), this phenomenon is due to the 
chordwise motion of either the transition point from 
laminar to turbulent boundary layer flow or a similar 
motion of the point of laminar boundary-layer separa- 
tion. 

By diminishing the theoretical damping by a rela- 
tively small amount, unstable oscillations also became 
possible if the flapping frequency was smaller than the 
pitching frequency and they showed good agreement 
with the observed oscillations. 

Finally, it may be added that unstable oscillations of 
rotor blades connected with vortex accumulation in the 
wake have also been observed by DuWaldt and Gates*® 
and Daughaday and Kline’ at Cornell Aeronautical 
Laboratory, but no theoretical explanation was added.t 
(2 THE VORTEX 


EQUATIONS FOR 


DISTRIBUTION 


} THe INTEGRAL 
We consider a particular blade, extending from —1 
to +1 along the X axis in the direction of the main flow. 
The the mid-chord of 
successive blades is 2h (see Fig. 2). 

On the blade, numbered n, which extends from —1 + 
2nh to 1 + 2nh, we denote the bound vortex strength 
by yn*(x, t) and the free vortex strength in the following 
If the vibrations are harmonic with 


distance between points 


wake by e,*(x, ¢). 
the frequency v we may put 
yn*(x, t) = Yyn(x)e”" (2.1) 
The relation between y*(x, /) and e*(x, t) is the law 
of conservation of vorticity 


¢ Since this paper was submitted, the authors have become 
aware of R. G. Loewy’s solution of the same problem by a 
different analysis (Journal of the Aeronautical Sciences, Vol. 24, 
No. 2, pp. 81-92, 144, February, 1957). Loewy’s results for the 
aerodynamic forces when h = 0 (no flow through the rotor plane) 
agree completely with ours. In the present paper results of 
flutter calculations with two degrees of freedom as well as a com- 
parison with experiment have been included 
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(Oy* Ot) + (de*/dt) = O (2.2) For harmonic vibrations, 


° * /(.. — Je 
In linearized theory, es, tf) = Ave 


there w = p Or: ws CO x Cc st; ; 0} 
(Oy* /dt) + Or(de* Ox) + (De* B!) = 0 (2.3) where w yi Qr and .1 is a complex constant. For eae] 
blade the value of .! is fixed by introducing the valy 


where & is the angular speed of the rotor and + the e, at the trailing edge 
distance of the section to the rotating axis. ; aici 
A; : p E, = t 
Since y* = 0 in the wake, we obtain 
y tvl —tw(Xx 1 2n/ 
eX, t) = €,°¢ 25 
(Oe* Ot) + Qr(dOe*/Ox) = O (2.4) 
The integral equation is obtained by calculating th, 
. | a -/Or ° ° 
or = Fo (x /Qr) | downwash due to the vortices and equalling this ex 
pression on each blade to the downwash, which is give 
by the motion w,(x.)e” 
»! kh = >—] 2} 2k wt 
| (Y; +- €; dé | ie 7) , . 
W,(x) = pe T ed EP - dé 
2m | J —1+42ki be 3 rr | J 142ki vr —& 


—|] + 2Pnh<v<l + 2nh (26 


If the mutual interference of the blades is neglected, the contributions of all blades with k # » and the wakes with 


k#n,k #n — 1 are omitted. This means that the influence of all vortices at a distance larger than 2/ of the mid- 
chord point is neglected. 
Then, the equation for x = 0 becomes (omitting the subscript 0) 

+] 2h —1 Jedd o—4 coal 
; fats | ete | e~ it 

+ tw - w(2h l % o 

WN) = = dé sia eC r dé + €_, ¢ : dé (2.3 
ZrJ-1 x —& 2a J x —& 27 J-2+1xn —€& 


For the other blades similar equations hold. 
If the contribution of vortices at distance larger than 
2h is neglected, we can write to the same degree of . ad 
ase and a bar denotes the complex conjugate. The solutior 
approximation . ‘ . 
of the integral equation 
*! 


(12x) | (vy + ©) /(x — &) Idé 
I 


7+ 
w(x) = (1/27) [(y + €)/(w — &) |dé 4 
I 


“ 
e 


Nou(x) — A-m(—xX) (28) 


iw2h iw v(x) = w(x) — Agqu(v) + A_yu(—~x) 


where No = @e’, Ay = €1€ (2.9) 


follows the same lines as the solution of the integral 


twee . t ie 
w(x) = flew — Ode = nthe iaysice theory 
J ies equation for a single airfoil (e.g., see reference 8). 


. ies Starting from the well-known solution of the integral 
wr wi Je /% (2 \ 4 : pena 
. 7 dg s) 10 equation for steady airfoil theory 
a e+] 
u(—x) = | [e'**/(—x — &) |dt (1/27) | I(y + €)/(x — &)|dé = v(x) (2.12 
J J —1 
—¢ | (edt /¢) we obtain, putting 
Jl x 
x =cosy, ££ = cos# 
2 ¢ [(" Wd)(1 + cos ddd 2 re i 2 . iia a(djdd 
yYy+e= tan = tan v(d)dde + sin ¢ - 
T S JO cos J’ — cos ¢ T Zz v0 T J0 cos? — cos ¢ 
») °n ») °n r « 
2 F u(d)dd Zz . ular — dV)dd " 
Ao sin ¢ + A, sin ¢ (2.13 
T J0 cos vd’ — cos ¢ T J0 cos vd? — cos ¢ 
From Eq. (2.3) it follows that on the airfoil 
lw(y + e) + (O€/Ox) = O (2.14 
or (Oc O¢) = tw(y + €) sing (2.10 
; : "? De iw (* l ; *™ o(9)(1 + cos ddd 
and, by integration, ¢(¢) = « + dy = e+ tan y sin ydy (2.16 
Jo oy J0 2 J0 cos 3 — cos wv 
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valuation of the integral gives 
| oe Lae 2iw [" sin (1/2)(8 — ¢) 
OF each e(g) = & + 2lw v(d)(1 + cos ddd + (3) sin d-In dd (2.17) 
7 J0 un J10 sin (1/2)(8 + o) 
© valu 
For ¢ = 7: e=e, WOT" = d_c™ 
and we find the relation 
25 ) ‘a 
; e’—he. = 2lw v(d)(1 + cos d)dd = 
ing the _ 
his e , a 9\ 7, ; ¥ ih 
“= 21w wid)(1 + cos ddd — 2iwdy u(d)(1 + cos d)dd + 2Ziwdr_, u(r — 8)(1 + cos d)dd (2.18) 
> 21Ver v7 /7 0 J 0 
Then we may write 
om ¢ ico 2iw [ sin (1/2)(8 — ¢) 
dv) = Xe 1 — + \ye*g + wd) sind In — di — 
~ r J0 sin (1/2)(8 + ¢) 
- Tw : oe sin (1/2)(8 — ¢) 2iw & sin (1/2)(8 o) 
2.0 No u(d) sin 3-In dy + Ay u(m — V) sin d-In | dy (2.19% 
a J0 sin (1/2)(0 + ¢) T J0 sin (1/2)(% + o&) 
S with : : 2 
d Reducing the integrals we find 
> mid- 
iw {" ; sin [(# — ¢) 2] Drv (" pwd) sin odd 2, (7 ula — BV) sin gdd 
et = wid) sin d-In — dy + — (2.20) 
Tr J0 sin [(o + oy) 2] r J0 cos ¢ — cosv Tr JA cos ¢ — cos J 
and subtracting this from Eq. (2.135) 
a8 
\o./ 
eo sin ¢ a sin (1/2)(8 — ¢) 
s(¢) = wd?) — wwsnd-In — dd + 
r J0 cos J) — cos ¢ sin (1/2)(8 + ¢)) | 
») o °n er 2 
tan | w(d)iddi — Ay ul(d)dd + A, u(r — d)dd (2.21) 
ution 1 2 LJo J0 J0 
lhe integrals over u(#) can be expressed in Hankel functions 
: sin ¢ os sin (1/2)(8 — ¢) 
vg) = wd?) — wsind-In -. dd + 
r J0 cos J) — cos ¢ sin (1/2)(08 + ¢) 
2.11 2 1 os ua (4) 4 
tan ¢ wid)dse — [AcZTo°? (w) + ATo (w) ] 2 (2.22) 
T 7 Jo | { 
oral 
The relation between the free vortex strengths at the leading and trailing edges becomes 
gral mn 
| wid)(1 + cos ddd = wr [Lh (w) + to? (w) | + rd [Zi (w) + Tp" (@) (2.23) 
J 
12 These formulas pass into the results for a single air 
foil, if A_) = O; hence, the results are the same, except and 
fora modification of the ?-function. 
Each airfoil gives a relation analogous to Eq. (2.25) yw) + Lo (ow) - 
and there are exactly NV relations between the NV con e ITO (wa) + Hy (w)|] = B 
stants \,, from which these constants can be deter : ‘ 
ined we get the system of equations from Eqs. (2.25) and 
(2.9). 
3) CALCULATION OF THE CIRCULATION FUNCTION 
13 a= ae, + Bex 
We consider a N bladed rotor, each blade has the an = ae + Be 
same amplitude distribution. Let the phase difference 
between two successive blades be 27(k/N), where k = 
14) | 91 N — 1, corresponding to N different modes of 
ie P = : N-1 
vibration. Then the ratio between the downwash at an = wén_1 + Ben 
5 ° ° Qritk/N) 
I) | successive blades is equal ton = e" : : : 
es oi cae : From these equations we can derive that 
Putting & for the first blade, 
6 & = 7 'a = 1 *e = i a 





e '*(4/n) i wid)(1 + cos 8)dd = a : F 
Jo and the solution is 
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/ ~ (2nik/N) — 
& = a/[a+8B ] 
Since for N blades 2Nh = 2rr we find cr 4 
2wh = (2arv/NQr) = (27/N)-(r/Q) 
a 
and ao = = ; - ; —(2ni/N) [k+(v/Q) 
(® + 1H) + (A + 1H )e . 
If the given velocity w(#) is expanded in a Fourier series 
wd) = a +2 Z. Ay, Cos nd 
n= 1 Fi 
we find 
4 (dy + a@)e i (2rik/N) tw2/ wid _=* } 
oo = = ; = ih ; > €1 = €yv-1 = @€ ~ - “" = €9€ ; 
7° + 21T,°?)) aa e ; 7, 4 UT) 
where A = (2/N)[k + (v/Q)] ,=o | 


The coefficient of tan (1 2)¢ in the distribution of the bound vortices then becomes 
The latter 


z j us a NM- { 
9 rac — ee [tH (w) + eo 4D (w) |¢ = snall unle 
T, 4 In the 
| H,®? + e—™ AH tHy® + eH yQ. The 
ys ly z , — @¢ 7 1 
"" h® + iH) +e ™(N® 4 iy) “1 EN® + iy) + em ™ (AN + iH") |] behavior | 
»Q = im 
Introducing the coefficient of ay in this expression as the P-function, ference” 
IT? + e™ WO unfavoral 
P(w, \) = ; a =F = aia N phenome! 
(HT, + 1H) + eo (© + 1H! 1 + 2[(Ao + e Hy) /(® + e7™™HM®) | men 
| 





1 + 7) [Jo(w) + Vo(w) tan (1/2)rA]/|Si(w) + Vil(w) tan (1/2)rd]} 


which is a circle in the complex P?-plane with center at 














P = 1/2 and passing through the points P = 0 and 
P= 1. Tand J 
- ° . ° = { ( ¢ 4 = ) es . 
We see that all expressions of the aerodynamic de P. [b+ O01/h) and Pa L + O(1/h) if inertia 
rivatives are the same as in the case of a single oscillat- Since we have neglected quantities of order 1/h, we see trifugal fe 
ing airfoil if we replace the ordinary P?-function by the that for the steady case the usual value P 1 is again | aerodyna: 
expression given here. obtained. mated by 
For a two-bladed rotor we consider symmetric and Furthermore. section 7 
antisymmetric oscillations. In this case, VN = 2 and, = (1/4) 
for symmetric oscillations, k = 0. . 
Then, R ) Ly = 1 
NUMBERS NEAR POINTS DENOTE jp 
\. = Q- — Op 06 | | | | [ ] 
Neg = V/X, W@W = V/ <M 
For antisymmetric oscillations, k = 1 and In = 1 
: 04 
Ae = 14+ (o/Q); w = v/Qr Lp ! 
In the first case, +0.2 where [ : 
P,(w, ») = pitching 
| Co io | ping ang 
-_ - : , ind the | 
1+7)[Jo+ Votan (rv 22)| [J + Vi tan (av 2) ]} : | 
, of the bl: 
(2.9) -0.2 
pore Moments 
and in the second case, ingles. 
-04 For tl 
) i] — 
P,(w, v) flutter cé 
l ave 
5 , - 06 4 4 1 | i { i } ' } have the 
1+2)[Jo+ Yocot (rv 20)]/[J; — Yi cot (rv 2Q)]] = | | | I= 0.54 
(3.3) 0 02 04 06 a8 i #£ C=, <« 
. : Fic. 3. The circulation function P,(\) for symmetrical oscilla- The ec: 
[If v/Q = 0, we obtain tions of a two bladed rotor (r = 15 
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i ROTOR - AX/S 
Fic. 4. The construction of the rotor system. 
P=0 if Si(w) + Vi (w) tan (mv/2Q) = 0 

or v/Q = 2n(n = 1, 2, 
» = () if Sy(w) — Yi(w) cot (rv/22) = 0 
or v/Q = 2n — l(m = 1,2 ) 


The latter approximations are allowed since w remains 
gnall unless v/Q becomes large. 

In the complex plane the vector ?(w) rotates with 
»Q2. The behavior is completely different from the 
behavior of the P-function of single airfoil. The case 
»Q = integer hence corresponds to a certain ‘‘inter- 
ference’ of the vibrations which, in 
ufavorable circumstances may be the cause of flutter 


and rotations, 


phenomena. 





Ve See 


agall 


}Z{1 — (Q2/y2)] + Ly 


cil — (2? 
{inertia with regard to these two axes. 
aerodynamic moments about the axes. 


section r = 


where / is the semichord and é/ the distance of the 
pitching axis aft of the quarter chord axis. The flap- 








illa- 





ping angle a is taken positive in downward direction 
and the pitching angle ¢ is positive if the trailing edge 
of the blade is more downward than the leading edge. 
Moments are taken positive in the same direction as the 
angles. 

For the statically balanced model for which the 
flutter calculations were made, the various parameters 
have the following values: / = 0.025 m, Rk = 0.5 m, 
I= 0.54 X 10-* kem:sec.”, J = 0.56 K 10—! kgm.sec.?, 
C = 0, ¢ = —0.1, while p = (1/8) kgm. —‘sec.?. 

The calculations were performed in such a way that 


HELICOPTER 


(ke/v*)y;a + Cll — (Q 
v)J + Lata + JJfl — (Q? 


/and J are the moments of inertia with regard to the flapping and the pitching axes, respectively. 
The terms containing 2 denote the additional stiffnesses due to the cen- 
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Fig. 3 shows the P,-function, Eq. (3.2), for w = 
(1/15)A (r = 15). 


(4) FLUTTER INVESTIGATIONS 
(4.1) Description of the System and Equations of Motion 


Fig. 4 shows the construction of the rotor system. 
Each of the two blades can rotate about a pitching axis 
AB and about a flapping axis CD. The flapping axis 
can be replaced by leaf springs of various strengths, 
providing stiffness against flapping, while the pitching 
motion takes place against the action of a torsion 
spring. The flapping axis CD is rigidly connected to 
the hub, which can itself rotate freely about an axis FF 
parallel to CE (the principle of the seesaw rotor). Asa 
result of the possible motion of the hub, the vibrations 
of the two blades are not independent, but can be 
separated into symmetric and antisymmetric vibra- 


tions. The axis FF is rigidly connected to the rotor 
axis. The blades are assumed to be infinitely rigid. 


The spanwise axis of each blade makes a small angle 
with the line EE. In this way the centrifugal force has 
a small chordwise component, which more or less com- 
pensates for the drag of the blades in rotation, thus 
reducing bending of the blades in the rotor plane. 

The equations of motion for flapping and pitching 
of the blade can be written in the following form: 


v)] + Lej¢e = 0 
(4.1) 


v”) | 7 Loe = (Rr rig = () 


C is the product 


trifugal forces. kg and ky are the mechanical stiffnesses against flapping and pitching. The terms Ly, etc., are the 
Their calculation involves a spanwise integration, which has been approxi- 
mated by the assumption that lift and pitching moment in each section are given by the value they would take in the 
(3/4)R when calculated with the theory of Sections (2) and (3). 


The blade extends between the sections 


=(1/4)Randr = R. Hence, the aerodynamic forces are given by 

In = Ly’ + tly” = mpl?[(3/4)R]*[1 (27P/w) | 

Li = Lye’ + ily” = mpl*[(3/4)R]?2) [(1/2) — e€] — (1 — ©€)(2iP/w) — (2P/w?) — (i/w)! 

In = Ly’ + ily” = mpl*[(3/4)R]?j [(1/2) — €] + €(2P/w) | aia 
Lon = Loo! + iL” rpl*{(3/4)R]} [(3/8) — « + &] + (e — &)(27P/w) + €(2P/w’) — | €)(1/w) | 


for a series of \-values, which determine at the same 


time w and P, the values of (vg/Q)? and (vr/Q)? for 
which harmonic oscillations are possible were calcu- 
/ / 
lated. Here vg = Vkp/I and vp = Vke/J are the 
natural frequencies for the nonrotating blade in flap 
ping and pitching. 
(4.2) Result of Calculations for a Statically Balanced 
Model 
The result of the calculations for symmetric oscilla- 
tions, using Eq. (3.2) for the P-function, is presented in 
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Fig. 5. For comparison, the results obtained by using 
quasi-steady aerodynamic forces (? = 1) and by using 
the aerodynamic forces referring to an airplane wing 
moving along a straight path are also presented. 

In Fig. 5 three regions of instability can be dis- 
tinguished. The region corresponding to the lowest 
values of \ (about 1.5) denotes a type of classical flutter, 
which in more or less good approximation is obtained 
by all theories, indicating that the wake is not of 
primary importance. 

The two other unstable regions are revealed by the 
new theory only. Here \ assumes values slightly above 
2 and 4+ so that the corresponding values of the P- 
function have small real parts and positive imaginary 
parts. The accumulation of vortices by the rotor 
rotating in its own wake is an essential condition for the 
existence of these unstable regions. The corresponding 
instability w.ll be called ‘‘wake flutter.” 

It is seen from Eq. (4.2) that for ? = 0, Li,” vanishes, 
which means physically the absence of aerodynamic 
damping for translational oscillations. Accordingly, in 
Fig. 5, the two unstable regions extend to (v7 Q)* = 

©, The corresponding values of (vz &)* are deter- 


mined from the condition 


T|1 _ (Q2? py”) | + Pg — (kp p”) = (4.3) 
which leads to 
(vp/Q)? = MCL + Jf apl2[(B/H)RE TI) -— 1 = 
1.024" — 1 


where A follows from the condition ? = 0 (yielding A = 
2.009 and 4.033). 

Similar calculations made for antisymmetric oscilla- 
tions would lead to regions of instability with \ slightly 
above 1, 3, etc. 

If the rotor speed is increased, a straight line toward 
the origin is covered by the relevant point in Fig. 5. 
For vz > vr the point would cross, successively, the 
various unstable regions beginning with the region of 
calculations showed 


largest \.. However, additional 


that the influence of structural damping is largest for 
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regions with large \. For the present system, a struc 
tural damping g = 0.01 made all unstable regions for 
which XQ is + or larger disappear. Experimentally, only 
the region with \ near 2 was observed. Experimental 
flutter was of a relatively mild type, which is in accord 
ance with the theoretical result that the size of the 


unstable region is sensitive to structural damping. 


(4.3) Influence of the Chordwise Position of the Inertia 
Axis 
Fig. 6 shows the result of the flutter calculations if 
C = +1 X 10 kgm:sec.?, 
inertia axis lying about | per cent of the chord aft of the 
Compared with the situation for the 


which agrees with an 


pitching axis. 
statically balanced blade, it is seen that the region of 
classical flutter has sharply increased, mainly in the 
direction of vyg~vr, which indicates the same type ol 
instability that would occur for a wing flying along a 
straight path. On the other hand, the regions of wake 
flutter are also enlarged, and those for \ = 2 and 4 link 
up continuously with the region of classical flutter. 


If the inertia axis is ahead of the pitching axis, wake 
flutter is the only instability which can arise. The 
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rr 


unstable regions are mainly limited to the part of the 
as is shown by Fig. 7 for C = 


fiagram where vg<r7, 
The wake flutter mode con- 


=6 X 10-* kem-sec.*. 
ssts of bending, to which some torsion is added for finite 

» 2)°. When the inertia coupling changes sign, the 
torsion is added with the same phase if vg> v7 is re 
placed by va< vr. 

Mathematically the imaginary part of the deter- 
minant of Eqs. (4.1) determines whether the unstable 
regions extend to plus or to minus infinity of (yp Q)? 
[his part is equal to 


Ata" Q?)(Q2 /p?) a IC{l — (92 py”) | + Lut La* 

tC[L — (Q?/v?)] + La’ j Li” | 
taking into account Eq. (4.3). Since Ly,” approaches 
zero from the negative side [see Eq. (4.2), Ke P>O], the 
sien of kp /Q? or (vp /Q)? will be negative if 


1C{1 — (Q?/0?)] + Li’ La” + 
Cll — (Q7/22)] + La’ fL” < 0 
According to Eq. (4.2), for P = C this becomes equal to 
Cll — (Q?/r?)| + La’ > 0 


Since in L.;’ only the aerodynamic inertia term 


mpl*|(3/4)R]7[(1/2) — e| 


is left, it means that the wake flutter regions extend to 


vp /Q)? = —o if the center of gravity of structural 4 


aerodynamic inertia is aft of the pitching axis. 
4.4) Influence of the Boundary Layer 

It was found experimentally that, for the blade with 
inertia axis (referred to structural inertia only) coincid- 


ing with the pitching axis, flutter could also occur if 
viz., in a symmetric mode with A~4 and in an 


Va< V7 
utisymmetric mode with A~3. These instabilities 
were not shown by the calculations (see Fig. 5). The 


mode consisted almost exclusively of pitching, suggest- 
ing one degree of freedom flutter. 

Previously’ it had been found, by measurements of 
aerodynamic forces on an oscillating airfoil placed in a 
normal two-dimensional flow, that the aerodynamic 
damping for pitching oscillations was considerably less 
than theory predicted if the Reynolds Number were 
near its critical value. In such a case there are transt- 
tion points at the upper and lower sides of the profile, 
Where the boundary layer changes from laminar to 


turbulent. The position of these transition points 


depends upon the angle of incidence and, hence, points 
motions during the 


Since the dis- 


there will perform chordwise 


oscillation of the profile (see Fig. 8). 
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placement thickness of the turbulent boundary layer is 
larger than that of the laminar layer and since the dis- 
placement thickness can effectively be added to the 
profile thickness, it will be clear that the originally 
symmetrical profile will become periodically slightly 
cambered. 

If the axis of rotation is before the transition point, 
the additional moment due to camber acts like a re 
straining spring. When a phase lag exists, the moment 
can be decomposed into a restraining component and 
an undamping component. The latter component is 
responsible for the overestimation of the damping in 
the theoretical coefficients. 

It has not been investigated whether transition of the 
boundary layer indeed occurs or that separation of the 
laminar boundary layer, which can act in a similar way 
as described above for the transition point, is the im 
The Reynolds Number in the 
reference section [ry = (3 4)R| is about 32,000, while the 
profile thickness is 12 per cent and the surface of the 


portant phenomenon. 


blade is rather rough. 
That the flutter phenomena for vz»: 
with boundary-layer 


vr are in fact 
closely connected effects is 
demonstrated convincingly by fitting a wire on each 
side of the profile at about 40 per cent of the chord. 
Flutter then disappears if yg<vr. The wire fixes the 
point of transition or of laminar separation 

One degree of freedom flutter becomes possible for the 
rotor if Lo.” mpl*|(8/4)R] [Eq. (4.2) | is increased by 3.5. 
This makes Ly.” positive for certain A-ranges. The 
value 3.5 is in good agreement with the differences 
between theory and experrments as measured in 
reference 5. 

Fig. 9 shows the result of flutter calculations near 
A = 4+ if Loy There 


now appears a flutter region for yy< vr also. 


is increased as mentioned above. 
There are 
two vertical asymptotes denoting purely one degree of 


freedom flutter. In this case, 


J{1 — (Q2/v*)] + Les’ + tLe” — (kr/s?) = 0 
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Fic. 9. Results of flutter calculations for statically balanced 
blade but with experimentally modified aerodynamic forces 


rpls (3/4)R] increased by 3.5 
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which yields \ = 3.58 and 3.97 as values for which Ly.” 


vanishes. The corresponding values of (v7/Q)? are 


(vp /Q)? = d*[1 + (Lm’/J)] — 1 


with 8.9 and 14.8 as solutions. 

Introducing structural damping by multiplying ky in 
(4.4) by 1 + zg, it follows that the damping for which 
harmonic oscillations are still possible is given by 


g = tan! (Ly»”/}J[1 — (1/d2)] + Lm’ } 


This value decreases with increasing J, which was in 
agreement with experiment since the symmetric mode 
for \ = 4 was replaced by an antisymmetric mode with 
\ = 3if J was increased above a certain limit. Appar- 
ently the structural damping then suppressed the \ = 4 
flutter but not the \ = 3 flutter. If J was increased 
still more, flutter became impossible for vg< vr. 


(4.5) Experiments 


Although a large number of experiments were per- 
formed, only a small series, which was specially suited 
for checking the foregoing results, will be described here. 
Experiments (a) to (e) were carried out with turbulence 
wires at about 40 per cent of the chord. 

(a) ve = 40.3 rad./sec., yp = 77.5 rad./sec., C = 0. 
No flutter. Agrees with Fig. 5. 

(b) vg = 40.3 rad./sec., vr 1s diminished to 26.2 rad. 
sec. by addition of weights near the flapping axis, C = 0. 
Hence, (vz /Q)? 


Flutter occurs for 2 = 18.3 rad./sec. = 
LS and (y7/Q2)* = 2.0, which is in the unstable region of 
Fig. 5. The experimental value of \ was 1.85. It is 
possibly due to a small rotation of the air with the rotor 
that this value is smaller than 2. 

(c) The blade was slightly overbalanced, C = —5 X 
10-* kgm./sec.?, vg = 36.6 rad./sec., yp = 25.1 rad. 
sec. No flutter. Agrees with Fig. 7. 

(d) The weights added in experiment 
C = —5 X 10* kgm.sec.?, ve = 36.6 
No flutter, although Fig. 7 


(b) were 


removed. 
rad./sec., vr = 66.0 rad. ‘sec. 
suggests that flutter should occur. Probably due to 


structural damping. 


(e) Overbalance increased to C = —-9 X 1076 
kgm.sec.”, vp = 35.6 rad./sec., vr = 61.8 rad./sec. 
Very mild flutter occurs for Q = 24.6 rad./sec., (vp /Q)? 
= 2.1 (v7/Q2)? = 6.3, ¥ = 1.94. This experiment 


shows that flutter becomes possible by overbalancing 
the blade, which is in agreement with Fig. 7. 

(f) Same experiment as (e) but without turbulence 
wires. Flutter occurs at 2 = 18.8 rad. ‘sec., (vp /Q2)? = 
3.6, (vr/Q)? = 10.8, X~4. 

(g) The overbalance of the blade was removed, but 
the natural frequencies were kept constant by adding 
weights. Flutter begins at Q = 30.3 rad./sec. with 
A~3 (antisymmetric). 

Since this is different from the result of (f), this 
means that overbalance favors this type of flutter, at 
least in the region vg< rr. 


(h) The additional weights also removed. 


were 
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Same situation as (a) but without turbulence wire 
Flutter at Q~21.0 rad./sec., with A = 4. 
firms the conclusion that a decrease of J makes flutter 
regions at larger \ also of importance. 


This con- 


(5) CONCLUSIONS 


The aerodynamic forces on oscillating helicopter 
blades without forward speed may be approximated by 
the usual formulas referring to a single airfoil in two- 
dimensional flow, provided a modified circulation func. 
tion is used. 

It has been shown both theoretically and experi- 
mentally that a rotor rotating in its own wake (during 
ground running or near autorotation) may experience 
unstable oscillations, which are generated by the regular 
vortex pattern of the wake (so-called wake flutter). 

If the inertia axis is aft of the pitching axis, wake 
flutter is mainly possible if the natural flapping fre 
quency is larger than the natural pitching frequency. 
For inertia axes ahead of the pitching axis the reverse 
holds. 

Wake flutter with potential theoretical aerodynamic 
forces can be eliminated by complete static balance of 
the blades. A backward position of the inertia axis is 
more critical than a forward position. 

At low Reynolds Numbers, boundary-layer effects 
may diminish the aerodynamic damping, thus promot- 
ing unstable oscillations even for a statically balanced 
blade. These can be avoided by fixing turbulence 
wires near the leading edges of the blades. 
In the case of 


However 
this is only of importance for models. 
full-scale helicopters, the Reynolds Number will be 
sufficiently large to prevent this type of instability. 
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Use of Stiffened Cylinders vs. Use of Unstiffened 
Cylinders 


Bertram Klein 

Chief of Stress, National Rocket Corporation 
Culver City, Calif. 

April 15, 1957 


I A RECENT PAPER,!' the authors presented an interesting com 
pilation of available test data for the compressive buckling of 
unstiffened circular cylinders with and without the presence of 
internal pressure. Certain curves were recommended for design 
Unfortunately, the inclusion of all test data lowers the previously 
proposed design curves significantly, even for the case of pres 
surized cylinders. Of course the pressurized cylinder is much 


more efficient than the unpressurized cylinder However, in 
many practical applications, ground handling loads are present 
while the cylinder is unpressurized. This condition may design 


the structure with consequent weight penalty Similarly, condi 


tions of intense acoustical loading or vibration, etc., may exist 

On the other hand, it is known? that the strength prediction of 
stiffened cylinders usually is much more reliable than that of 
unstiffened ones and that the strength-weight ratio of efficiently 
stiffened shells exceeds that of unstiffened shells. It is believed, 
therefore, that weight saving is possible in many practical cases 
if stiffened cylinders are used in place of unstiffened although at 
an inerease in cost and complexity of fabrication. The actual 
form of the stiffening may vary from application to application 
It appears that stubby integral stiffening, for example, may prove 
useful when thermal effects are present. Sandwich construction 
of various types may be considered to be a form of stiffening but 
with less strength reliability than integral stiffening at least at 
the present stage of development 

In conclusion, it is advocated that more emphasis be placed 
on the development and strength prediction of stiffened cylinders 


rather than unstiffened 
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On the Velocity Profile for Turbulent Flow 
Near a Smooth Wall 


John W. Miles 

Professor of Engineering, University of California, Los Angeles 
Calif 

April 29, 1957 


poe DrieEst! RECENTLY has derived a continuous expression 
for the mean velocity profile in turbulent flow near a wall 
on the hypothesis that Prandtl’s mixing length should vanish 
like 1 exp(—y/A) as the wall (y = 0) is approached. The 
result 
in the boundary layer if the parameter A is chosen to produce 
an asymptotic fit over the logarithmic portion of the profile 

The writer arrived at rather similar results by m«cdifying 


closely approximates experimental data everywhere 


the mixing length should be pro 
to the hypothesis that it 


Prandtl’s original hypothesis 
portional to the distance from the wall 
should be proportional to the distance from a sublayer of com 
pletely laminar flow. This is perhaps, on physical grounds, less 
satisfying than the mechanism of damping proposed by Van 
Driest, but it leads to a somewhat simpler result that fits the data 
equally well 
We start from Prandtl’s expression 


t/p = v(On/Oy) + 1(On/dy)? (1) 


for the total shearing stress 7, where p, v, 4, and / denote the 
density, coefficient of kinematic viscosity, mean velocity, and 
mixing length. Weassume that 


l= K(y— yw), y= ! 
= () Of yS yo) 


’ 


‘e 


(2) 
where A is von Karman’s universal mixing constant, and yo de- 
notes the edge of the lamittar sublayer. Introducing the dimen- 


sionless quantities 








ue =(r/p) "7a and ye = (7/p)!*y/p (3) 
we obtain the solution 
ye = yor + (1/2K) sinh a (4 
use = yor + (1/K) [a — tanh (a@/2) (5) 
Oux Oy = sech? (a/Z) (6) 
16 See > 
Ux = You \/ | m 
° | 
|2+——-+ . = ——- : t 
4 | | — VAN DRIEST, 
a) a | | ‘. K=0.4, Ay=26 
amen MILES, 
| | K=0.4, C=5.24 
| SS ae | | +. LAUFER DATA 4 
Re 
° 50,000 
| * 500,000 } 
O 10 20 30 40 50 60 70 
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with @ as a parametric coordinate for y > yo; for y < yy we hav 
the laminar solution wx = yx. We remark that (Ou /dy, 
is the ratio of the viscous shearing stress to the total shearing 
stress 

The parameters A and yo may be evaluated by comparing th 
asymptotic representation of Eqs. (4) and (5)— viz., 


ux = (1/AK)in ye + yor + [(In 4A L)/A " 


to an experimentally determined profile in the form 
we = (1/K) In ye + C x 


We adopt the value A = 0.4; comparing Eq. (7) to Eq. (8) th 
vields 


yor = C + 1.825 q 

The profile given by Eqs. (4) and (5) for C = 5.24 is plotted jy 
Fig. 1. It lies slightly above van Driest’s curve and fits Laufer 
data for Re = 50,000 (as presented in van Driest’s paper 
somewhat more closely; on the other hand, it gives a rather less 
satisfactory fit to the Re = 500,000 data than does van Driest’s 
result 

We conclude by noting that, on the basis of the above model 
larger values of C imply thicker laminar sublayers. This is of 
some interest in connection with measurements made. over 
smooth water; thus, Roll’s measurements? for (7/p)!/2 < 10 em 
sec. lead to values of C considerably larger than 5.5. Prandtl! 
inferred from this that the flow was laminar throughout the first 
2m. of the boundary layer, but this is difficult to accept, and the 
explanation of an increase (greater than would have been anti 
cipated on the basis of smooth wall data) in the thickness of the 


laminar sublayer appears to be more attractive 


ADDENDUM 


In reference to a recent note by Elrod,* it perhaps should be 
emphasized that the foregoing considerations are not intended 
to predict the manner in which the turbulent shear stress falls 
off very close i.e., ¥ < Yo to the wall 
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On a Minimum-Time Flight Path of a Jet 
Aircraft 


John Carstoiu 
Associate Scientist, Scientific Research Staff, 
Republic Aviation Corporation, Farmingdale, N.Y 


April 25, 1957 


SECTION (1) 


— PROBLEM of minimum time to climb of a jet aireraft has 
lately received widespread attention. Although theoretical 
and experimental data on this problem have appeared in several 
publications and reports, f the subject is still apparently of a con 


troversial nature and has led in the literature to some misin 


+ We are indebted to Dr. van Driest for the reproduction of Fig. 2 from 
his paper 

tt Analysis compiled under the direction of Dr. Theodore Theodorsen at 
Republic Aviation Corporation 

t Only a few references are quoted here 
references given by Behrbohm! and Cicala and Miele.? 
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SECTION (2) 


The problem can be approached by considering the motion, in 
vertical plane, of the mass center of the aircraft, having a mass 
jal to that of the aircraft and acted on by a force equal to the 
ctor sum of the external forces. These forces are the in- 
‘antaneous weight of aircraft W, the net thrust 7, and the aero 
namic forces—drag D and lift L 

Consider the moving frame, with respect to 
h) supposedly fixed, defined by the positive tangent 


a rectangular 
rame (¥X, 
dnormal to the path of the mass center 

By assuming that the thrust acts along the tangent to the path, 
emotion of the mass center is described by the equations* 


(W/g)V = T — D — Wsin y) 
(W/g)Vy = L W cos 4 \ 


S 


(2.1) 


ere V is the speed, g the acceleration of gravity, 7 the positive 
dination of the path to the horizontal, and the dot indicates 
rivatives with respect to the time 
In addition, one has to consider the rate of change of weight due 
fuel consumption. This can be described by the equation 
W = —zl (2.2 
ing the specific fuel consumption 
The path is completely determined by the integration of the 
juations 
t = V cos y/ 92 
h = Vsin y\ 


[he time to climb is given by 


° 
. . , A ” 9 
time to climb = f (ds/V), ds = YUdx? + dh? (2.4) 

along 
path 
One has that 
D = (1/2)pSV2Cp, L = (1/2)pSV2C,_ (2.5) 
yhere the nondimensional coefficients Cp and Cy, are functions 

{Mach Number M and the angle of attack a@ (positive ineli 

ition of datum wing chord to direction of flight). Furthermore, 

ne has the polar’s equation 
Cp = Cp, + Co, (2.6) 
vhere the induced drag coefficient Cp; admits, generally, a rep- 


resentation of the form 


Cp = K(M, Cr) (2.7) 
One can similarly write 
T = (1/2)pSV2C7 (2.8) 


where Cr is a function of Mach Number and altitude only for a 
fixed power setting 


The specific fuel consumption Mach 


w is a function of 
Number and altitude 

Under these conditions, the differential Eqs. (2.1) and (2.2) 
in, theoretically, at least, determine V, y, and W as functions 
f time if the values of these variables are known at some time 4 
ind the angle of attack is specified. Eqs. (2.3) and (2.4) will 
then, respectively, determine the path and the time to climb 





SECTION (3) 


The problem is to minimize the value of the integral (2.4) under 
ettain boundary conditions. As the time ¢ is precisely the vari- 
ible to be minimized, one has to choose some other variable as 
inindependent one. Because all experimental data are functions 
f altitude (and speed), it is convenient to take f as the inde 
pendent variable. With h as the independent variable, Eqs 


2.1) and (2.2) can be written in a suitable form for numerical 








* The equation which states the equilibrium of pitching moments on the 
tircraft has no immediate relation to the performance problem, See, for 


nstance, reference 6 
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integration as follows: 


V’ = N(Cr — Cp) (g/V) ¢ 31 
y’ = (N/V)CL — (g/V?2) cot 75 rr 
WwW’ = (NW'/g)uC7 (3.2) 


where the prime indicates derivatives with respect to 4, and 
N = gpSV/2W sin 4 
is connected? with the unit of aerodynamic time ¢ by the relation 
2isiny = N 
Eqs. (2.3) reduce to 


x’ = cot 7 (3.3) 


The time to climb between two altitudes hy and h, is now given 


2h, 
} (dh/V sin > (3.4) 
h 


by 


The problem of minimum time to climb can be formulated 
as follows: Select from the solutions V(h), y(h), Wh), x(h), and 
a(h) of Eqs. (3.1)-(3.3) together with (2.6) those solutions which 
have initial given values Vo, yo, Wo, Xo, ho and given end values V,, 
Yer Xe, h, and which yield a minimum value for the integral (3.4) 


According to the appropriate formalism of the calculus of vari 
ations,’ the integrand in Eq. (3.4) is to be replaced by 


F = (1/V sin y) + Ay[V’ N(C7 Cn) + (g/V)| 4 
ryly’ (N/V)C. + (¢/V?) cot y] + An 
] 


[W’ + (NW/g)uCr|] + A,[x’ cot y]| (38.5) 


where Ay, Ay, Aw, and A; functions of altitude, are the Lagrange 


] 


multipliers corresponding to Eqs. (3.1)-(3.3). 
The Euler-Lagrange equations give 


\'y = Av IC N/V) (C7 Cp) + N(O/OV ) (C7 Cp) + 
(g/V?)] — (2d yg/V*) cot y + Aw(NW/g) X 
[(1/V )uCr + (0/0V) (nCr)]) (1/V? sin 
A’ = Av(Cr Cp) N cot y + (A,/V) X (9.6) 
[NC, cot y -— (2/V sin? y)| Aw X 
(NW/g)uCr cot y + (A,/sin? y) — (cot 7/V sin y) 
Nw = (N/W) [Av(Cr — Co) + (A,/V)Cz) 
A, = constant, OCp/OC, = (1 VQ, Av) 


which have to be considered together with the boundary condi 
tion? 


(3.7) 


[bvAV + A,Ay + AwAW +4 2, Ax] fe = 0 


If one does not specify the distance x to be covered, then one 
necessarily has 
‘4, =0 
On the other hand, the weight W cannot be specified at h = h,; 
it follows that 
Aw(h,) = O (3.9) 
Under these conditions, the problem has a solution. <A step 
by-step numerical integration is indicated, to be handled by a 


high-speed digital computer 
SECTION (4) 


If the operation involved does not require assigned values for 

V and y at the end of the climb (unrestricted or free climb), then 
one has the following conditions: 

Av(h.) = O (3.10) 


y(h.) = 0 (3.11) 


The angle of attack a,, at the end of the climb is determined 
by applying l’'Hospital’s rule to the right side of the last Eq 
(3.6). One has immediately, by virtue of the first two Eqs. (3.6 
and conditions (3.9)-(3.11), that 


= cot 


(OCp/OCL) 
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The Secondary Flow Behind a Cascade 


S. Soundranayagam 

Mechanical Engineering Laboratory, The English Electric 
Company Ltd., Whetstone, England 

May 1, 1957 


r REFERENCE | the vortex filament concept is applied to de 
rive expressions for the secondary circulation downstream of 
an impulse cascade. This concept can be extended to the case 
of diffusing or accelerating cascades, and it can be shown that 
the resultant secondary circulation is zero for an impulse cascade 

Let ea be a typical vortex filament in the entry flow, which is 
caused by velocity variation in spanwise direction such as a 
boundary layer. Let Ad be the position of this filament after a 
time ¢ has elapsed. Consider the motion of a particle at a and 


ate. 
t = (ab/U;) + SL qu) + (ed/U2) = : 
(ef/U,) + S (Lu qu) + (gh/ U2) 

where 

L,; = distance along lower side of blade 

L, = distance along upper side of blade 

g. = velocity along lower side of blade 

gu = velocity along upper side of blade 

U, = inlet velocity to cascade 

U. = velocity downstream of cascade 


i(ab — ef) U;} 4 i Sh qi) - Fu. qu)} = 
(gh — cd)/Us = dk/U, 


s(sin ay/U,) + Sgt gi) — Fh qu} = 


dk/U (cos a;/cos as) 


where s is the pitch, since by continuity U; cos ay = Ls cos az. 
dk = s(sin a;/U;) + iS (Li qi) - S (Lu qu) x 
U,(cos a;/cos a») (1) 


Consider two particles at e, one just above and one just below 
the stagnation streamline. The particle just above travels to h 
via the upperside of the blade while the one below travels to j via 
the underside of the blade. This results in the vortex filament 
being stretched. One can think of the vortex filament as being 
stretched directly from / to 7 or being stretched from h over the 
upperside of the blade to f and round the underside to ¢g and j 


t = (ef/Ui1) + J (Lu/qu) + (gh/U2) = 
(ef /Ui) + J(Li/n) + (¢j/U2) 
(gh — gj)/U2. = SS (Li qi) — SF qu) 
jh = iS (Li qi) — S Lu/qu)} (U1 cos a1/cos as) 








Consider an elemental circuit moving with the fluid enclosing 
an area dS which lies in the plane determined by a streamline ar 
the z direction (spanwise direction). Let its trace be as show 
in Fig. 1. The subscripts 1 and 2 refer to conditions before an 
after the cascade. The circulation dK associated with this cir 


cuit is 


dK, = wdS,; where w, = dU,/dz 
The elemental area downstream of the circuit will be 
dS» = dS\(cos ay/COS ae) 


Let the resultant vorticity downstream be w: and let it hav 
components @2, and a2. normal and parallel to the streamlines 


respectively. 





wendSe = wdS, 
by Kelvin’s theorem. Therefore, 
Wen = w1 (COS a2/COS a1), wes = way, COLO 
where @ is the angle shown in Fig. 1. 

wes represents a streamwise vorticity which is distributed right 
across the exit passage. The secondary circulation associated 
with this distributed vorticity is 

I; = id- (unit distance in z direction) (average value of w 

= 1d-wo, (average value of cot @ over the passage) 

If the curve of the vortex filament dh is of moderate camber 
the average value of cot @ can be shown to be cot ¢, where ¢ is 
the angle between the streamlines and a straight line joining dh 

I; = td-an- cot d 
= 1d-w(cos a2/cos a)-(ih/id) 


= (COS ae/COS ay )ih 


th = dk s sin a 


I; = o(cos a»/cos a) [dk 5 SIN ae 
which gives 
Ir, = (dU,/dz) |\(s/2) }(sin 2a; sin 2a2)/cos a} + 
tS (Li/q) — J (Lu/qud} 3 
The direction of T; is that of a right-handed screw facing 


upstream 

The component of the downstream vortex filament jh repr 
sents a vortex sheet springing from the trailing edge. Consider 
the particles at # and j to be separated by a very small distance 
dn normal to the streamlines. The filament is considered to be 
stretched from 7 to h. 

As before, 


wes = Wey Cot ¥ 
= won(jh/dn) (since dn is small) 
= wi (COS ax/cos a) (jh/dn) 


The circulation in the vortex sheet of unit distance in the : 
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fjrection is 
Is = wi(cos a2/cos a) (jh/dn)dn 


COS Qa ) jh 


= wWCOS a 


which gives 
(dU,/dz) {SS (Li qu) SL qu) Uy (4) 


The circulation of this is a right-handed screw facing down 


> = 


stream 
The vorticity due to the shed circulation from the blades is 
iso present downstream The circulation associated with this is 


Ir; = (dU,/dz)s cos aj|tan ay tan a»| (5) 
The direction of T’; is the same as Ty 
The total streamwise circulation is therefore 
r=T) rr — V3 
T : ‘ : ‘ ‘ t 
= s(dU/dz) [}(sin 2a sin 2a2)/2 cos ay} 
cos a;(tan ay tan az) (6) 
when a} = —ae as in an impulse cascade 
r=0 (7) 


The above method may be applied to an accelerating cascade 
when the following results are obtained. These results are the 


same as those obtained by Hawthorne? with a different method 


r, = (dU; /dz) (s 2) {(sin 2ax — sin 2a:)/cos a} 4 

tS (Li/qr) — J La/qu)i 
r = (dU,/ds) | f(Li/q) — J (Lu/qu)| 
Tr; = (dU,/dz)s cos a,(tan ae tan a) 
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On Ménard Inserts in Supersonic Nozzles} 
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INTRODUCTION 
Intentional changes of the design Mach Number in a nozzle 
with fixed contours by means of movable sidewall inserts were 
first reported by Ménard.' His work formed the starting point 





for studies in the Aerodynamics Laboratory of The University 





‘ing 


re 
der 


ice 


of Michigan into the potential usefulness of such inserts at low 
and high Mach Numbers. The study was motivated by the 
possibility of using Ménard inserts as a simple and extremely 
economical technique of changing the Mach Number of a set of 


fixed nozzle blocks 
DISCUSSION 


Consider 
uniform exit flow of Mach Number 9; add to the plane walls of 
this nozzle properly shaped and located half bodies such that a 
new uniform flow of Mach Number 14, = JJ) + A.V results 
The determination of the contours of these inserts and of their 


a given two-dimensional symmetric nozzle with a 


location relative to the basic nozzle cannot be carried out by any 
of the usual design procedures for two-dimensional or axisym 


Only the simple one-dimensional area ratio 


metric nozzles 
t This work was part of the Bumblebee research program at the Engineer 
ing Research Institute and was supported by the Bureau of Ordnance under 
Contract NOrd 16595 

* Research Associate 
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*** Associate Professor of Aeronautical Engineering 
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theory or the cumbersome technique of characteristics for gen 
eral three-dimensional supersonic flows appear applicable. The 


specific design procedure employed is given in references 2-5 


RESULTS 
Low Mach Number Inserts 
The technique described above was used to design a set of 
inserts for the 19 in. X wind tunnel of the Ordnance 
Aerophysies Laboratory so that their Mach 1.5 nozzle would give 
Mach 1.6 flow of acceptable uniformity 
Inserts of the recommended design were manufactured and 
The results 


27.5 in 


tested by the Ordnance Aerophysics Laboratory. 
were highly satisfactory and are reported in reference 2. Briefly, 
the desired Mach Number of 1.6 was achieved with the inserts in 
design position, and the flow uniformity was no worse than that 
of the basic Mach 1.5 nozzle. Off-design positions of the inserts, 
6 in, upstream and downstream of the design, resulted in usable 
flows of Mach 1.57 and 1.63, respectively, with the most uni 
form flow being obtained in the latter position (Fig. 1). The 
cost of these inserts and their installation amounted to approxi 
mately 13 per cent of the expected cost of a new set of Mach 1.63 
nozzle blocks. 
High Mach Number Inserts 

After this success at low Mach Numbers, the program of d« 
signing inserts for the OAL facility became quite ambitious. The 
first high Mach Number inserts were to boost the OAL Mach 
2.23 nozzle to Mach 3.0. It was found that no inserts could be 
designed which would have a reasonable chance of giving ac 
ceptable flow at the Mach 3.0 level.® 

Next, a series of inserts (Fig. 2) were studied to bring the OAL 
Mach 2.75 nozzle to a Mach Number of 3.0. Exploratory tests 
with models of the proposed inserts were conducted in the 
Mach 2.84 nozzle of the 8 in. X wind tunnel of the Uni 
The first series of tests were made of so-called “short”’ 


13 in 
versity. 
inserts, shown in the upper portion of Fig. 2. The negative 
results of these tests are given in reference 4. The unsatisfactory 
flow of ali three short inserts was primarily due to a strong lateral 
shock believed to originate in local regions of rapid expansion and 
cross flow just downstream of the effective throat of the nozzle- 
insert combination. 

In order to alleviate the local expansion, a long cylindrical 
insert was tested® as shown in the lower half of Fig. 2. The main 
result is given in Fig. 3 and shows that a flow of reasonable uni 
formity was obtained, and that the Mach 3.0 level of the basic 


nozzle was increased by an increment of 0.20. The lateral shocks 
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Mach Number (Univ. of Mich. Mach 2.84 nozzle). 


mentioned before were now sufficiently weak and were quickly 
attenuated. 

From these tests one can predict a usable Mach Number 3.0 
flow for the OAL Mach 2.75 nozzle with cylindrical inserts 
However, in case one wants to terminate these cylindrical inserts 
in a smooth trailing point upstream of the test section, additional 
model tests are strongly recommended to assure that the flow 
uniformity remains acceptable. 

CONCLUSIONS 

The work with Ménard inserts to date has indicated the fol 
lowing: 

(1) Ménard inserts are an effective and economical means of 
increasing the Mach Number of an existing two-dimensional 
nozzle. 

(2) At low Mach Numbers (around 1.5), increments of up to 
10 per cent of the basic Mach Number can be achieved without 
loss in flow uniformity. 

(3) The design of inserts for the low Mach Number range 
can be carried out by a modified area-ratio theory. In most 
cases no model tests will be needed. 

(4) At high Mach Numbers (above 2) increments in Mach 
number due to inserts will be less than 10 per cent and some 
deterioration of the uniformity of the flow must be anticipated. 

(5) The design of inserts for the high Mach Number regime 
should be substantiated by model tests, especially in the case of 
large increment inserts in short nozzles. ; 
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Physical Mechanism of Vortex-Ring Cascade 


Norbert J. Dabrowski 
Research and Development Department, Electric Boat Division, 
General Dynamics Corporation, Groton, Conn. 


May 10, 1957 


y A RECENT LETTER,' C. S. Hsu reports the observation of a 
vortex-ring cascade phenomenon. The writer suggests that 
the phenomenon may be explained by the following physical 
mechanism. 

Consider a dye drop entering a quiet body of fluid and forming 
a vortex ring. Assume the dye to be denser than the surround 
ing fluid. The vortex ring would be a perfect torus cnly if the 


dye drop were perfectly symmetric and the fluid absolutely 
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tortion is upward [Fig. 1(b)]? 


SEPTEMBER, 1957 


quiet. Such a ring would be in unstable equilibrium and so woy| 
not initiate a cascade as long as the motion and conditions re 
mained ideal. Practically speaking, however, the ring will hay, 
small distortions. For a downward dip [Fig. 1(a)], the pressyy, 
at the level of the undistorted part of the ring is the same in th, 
fluid and within the dye filament, but along the dip, this equalit 
does not hold. Since the filament is continuous, the pressuyr 
within it at the bottom of the dip is wah (dye density X head 
just as it would be in a mechanical tube. The pressure in th 
fluid at the h level is w,h (fluid density X head), which is less 
than that in the dip since wy < wy. Therefore, the dye in th 
lower part of the dip tends to push away the fluid immediate} 
surrounding it. Since there is no restraining tube (other than th, 
effect of surface tension), the dye will flow from ring level dow; 
into the dip thus expanding the dip [Fig. 2(b, c)]. 

Meanwhile, the entire vortex ring is moving downward. Thy 
drag force created by this motion causes a flattening at the bot 
tom of the dip such that the dip takes on an approximatel 
ellipsoidal shape [Fig. 2(c)].. The asymmetry, which is bound t 
be present in the ellipsoid, will allow the entering filaments t 
take on a configuration that produces a rotational flew about 
vertical axis [shown by the arrows in Fig. 2(c)]. By centrifuga 
action, the ellipsoid changes to a torus. The central part of th 
ellipsoid gets thinner, becomes a membrane, and breaks [Fig 
2(d, e)}. The vortex ring, which has been supplying the dye for 
the small pro forma rings, has emptied; the near-vertical filaments 
break [Fig. 2(e)], leaving a torus. While the rotational flow has 
been growing, the pro forma ring continues to descend. The 
relative upward velocity of the fluid past the pro forma ring 
[Fig. 2(b, d, e)] imparts a vortex motion so that after the mem 
brane has broken a vortex ring results [Fig. 2(f)]. 

Returning to the original vortex ring, what happens if a dis 
At the top of the distortion, th 
At the level of the 


undistorted ring, the pressure within the filament is more than 


pressures in the fluid and the dye are equal. 


that in the fluid because the dye density is greater than the fluid 
density. Therefore, the dye at the ring level tends to push away 
the surrounding fluid, thus permitting the dye in the distorted 
portion to flow down to ring level and so reduce the amplitude of 
the distortion. It can be seen that the ring is stable to upward 
disturbances but will produce rings if given a downward dis 


tortion. 
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b. UPWARD DISTORTION 


Pressure conditions for distorted vortex ring when dye 
is denser than surrounding fluid. 
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Formation of vortex ring from a downward distortion 
when dye is denser than surrounding fluid 


Fic. 2 





For a dye ring that is lighter than the fluid, similar reasoning 
shows the ring to be stable to downward disturbances but prone 
to produce rings when distorted upwards. Fig. 1(a, b) applies 


to this case if the subscripts are reversed. 
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Flow Over Blunt-Nosed Bodies 


Robert W. Truitt* 

Virginia Polytechnic Institute, Blacksburg, Virginia 
May 14, 1957 
NREFERENCE 1, experimental values of the dimensionless 
sagnation velocity gradient for a hemisphere-cylinder were 
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presented for a rather complete range of subsonic and super 
sonic Mach Numbers. Korobkin and Gruenewald! found that, 
due to the influence of the cylindrical afterbody, their subsonic 
experimental data agreed more with potential flow around a 
Blasius half-body than with flow around a sphere. As stated in 
reference 1, the half-body was selected for comparison at J, = 
0) since, to their knowledge, the problem of the hemisphere 
cylinder in potential flow had never been solved mathematically. 
It is interesting to note, however, that the hemisphere-cylinder 
and several other related bodies in potential flow have been 
solved in references 2, 3 and 4, by Dr. J. F. Vandrey, now with 
the Martin Company. 
tential flow problem of a cylinder with hemispherical, ogival and 


In reference 2, Vandrey solved the po 
flat rounded-off heads by the vortex-sheet method. Reference 
3 deals with the transverse flow past a hemisphere-cylinder as 
generated by source distribution on its surface. In _ reference 
4, Vandrey considers the case of a finite body with a hemi 
spherical head, a cylindrical central part, and a boat tail for 
longitudinal and transverse motion and for rotation about an 
arbitrary axis. Since reference 4 also uses the source distribu 
tion method, two independent calculations for the hemisphere 
cylinder are available. The results were found to be in excellent 


agreement with experiment. Vandrey’s results yield, for the 
hemisphere-cylinder, (8D/U.) = 2.8 and for the two-dimensional 
It is 


clear from Vandrey’s solutions that the influence of the afterbody 


cylindrical nose with straight afterbody, (8D/U’.) = 3.65 
is significant. It is interesting to note that (6D/U’ = 8/3 for 
the Blasius half-body! is in excellent agreement with Vandrey’s 


result for the hemisphere-cylinder 
PERFECTLY FLAT-NOSED BODIES 


Evidence of afterbody influence on the stagnation velocity 
gradients was also found for perfectly flat-nosed bodies in refer 
ence 5, based on the method of reference 6 
The condition for modified Newtonian-type pressure distribu 
tion to exist on the surface of cylinders and spheres was derived 
in reference 5 in a somewhat more satisfactory manner than was 
originally presented in reference 6. Using Eq. (11) of reference 
6, the ratio of local to stagnation pressure coefficient may be 
written (in notation of reference 6) 
Cye/Comsinue = 1 — (Ki /d)sin? 0 (1) 


De 
where 


h=1 + 32[((M,/M,)? 11/1 + yM,?)} (2) 


and K; is, for the sphere, 
K; = K; = (1 + } [2 + (b/a)*)/2[(b/a)* 1}{)? (3) 
and for the cylinder 


\(b/a)? 1}})? (4) 


K; = Ko = (1 + } [(b/a)? + 1) 


The condition for Newtonian-type pressure distribution to 


exist (see Eq. (1)) is simply 
K,; A= 1 (5) 


The condition shown in Eq. (5) not only yields for a given Mach 
Number M7 « 
b/a (for cylinder or sphere), but also predicts the minimum 
Mach Number for which the modified Newtonian pressure dis 
The calculated minimum Mach Numbers were 


the corresponding shock detachment parameter 


tribution exists 
approximately 1.38 for the sphere and 1.9 for the cylinder 

Corresponding results for the normal flat plate were found in 
reference 5 by methods of conformal mapping using the previous 
solution for the cylinder.6 By employing superposition, ap 
proximate solutions for two other perfectly flat-nosed bodies 
were found, namely, the blunt-nosed body with straight after 
body and the flat plate with free streamlines.‘ 

The results of the calculations for the dimensionless stagnation 
velocity gradient for the several blunt-nosed body types cited in 


the above references are presented in Fig. 1 
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On the Normal-Force Distribution of a Cone- 
Frustum Body at Supersonic Speeds 
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May 15, 1957 


HE SELECTION of a missile configuration depends upon a 
es of considerations, one of the most important being 
the aerodynamic characteristics. Conical-nosed bodies, par- 
ticularly conical-nosed bodies followed by frustum afterbodies, 
are useful missile configurations. In order to reduce the number 
of wind tunnel tests, an accurate theoretical method for pre- 


dicting these characteristics is very desirable. 
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The well-known perturbation theories can be used for the cone 
frustum body, but they are laborious to apply without the aj 
of high-speed computers. Van Dyke’s second-order and hybric 
theories are generally more accurate but they, like the first order 
theories, cannot be applied at combinations of Mach Numbe, 
and body shape for which the Mach cone is less than the tangent 
cone at the vertex. In reference 1, a correction is applied to tly 
hybrid theory to give the exact axial- and cross-flow values o; 
the initial nose cone. 

Recently Syvertson and Dennis? developed a_ second-order 
shock-expansion gives satisfacter 
results for bodies similar to the cone-frustum and for whic} 


theory which apparently 


computation is relatively easy. The theory can be used in regions 
where perturbation methods are not applicable 

While all the theories mentioned have given good results under 
certain conditions, no report of their specific application to 
cone-frustum body has been found. Thus, this present work was 
initiated and has resulted in the determination of the preferred 
theory from the standpoint of accuracy and ease of applicability 
This note presents experimental data and results of theoretica 
calculations for a cone-frustum body for a Mach Number rang 
of 1.54 to 4.38. 

The experimental model consisted of a 45° cone followed by 
10° frustum afterbody. Fifteen pressure orifices were mounted 
along the length of the model perpendicular to the exterior sur 


face. Pressure measurements were made at nine roll angles and 
for 0°, +2°, 4°, and 6° angle of attack. Local normal fore: 


coefficients were determined by integrating the pressure dis 
tribution at each station. These values were then used to deter 
mine total normal force and moment coefficients. Center of 
pressure was computed by the relation Y,,/D = Cye/Cva, wher 
D is the base diameter and Cyq and Cyq are the initial moment 
curve slope and initial normal force curve slope, respectively 
Fig. 1 presents local normal force coefficient distributions ob 
tained by theoretical methods and by experiment. These 
distributions are typical for the Mach Number range of 1.54 
to 4.38 


hybrid theories were carried out on an IBM computing machine 
: § 


The integrations for first-order, hybrid, and modified 


since the hand computational procedure is extremely laborious 


and time consuming 























MODIFIED HYBRID a - 
COINCIDENT WITH M*2.44 Q=2 | 
SECOND — ORDER © O EXPERIMENTAL DATA 
Sn  SHOCK-EXPANSION —> —— MODIFIED HYBRID 
te Xx ---- HYBRID 
i —-— FIRST-ORDER 
0s fi B —- MODIFIED FIRST-ORDER | 
~ fy —— SECOND-ORDER SHOCK-EXPANSION | 
fA. : J 
Ce ee ee 
el 
x 
y 
20 ie 

















Fig. 1. Theoretical and experimental normal-force distributions 
on a cone-frustum body at Wo = 2.44 
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Fig. 2. Initial normal-force slope versus Mach Number 
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ment 
y It can be seen from Fig. 1 that first-order theory vields inac 
$e urate results, underestimating the cone loading while overesti 
boars mating the loading on the frustum. Perhaps even more im 
sii portant is the fact that the shape of the normal force distribution 
lified ishighly inaccurate causing large errors in estimation of the loca 
hine tion of the center of pressure. Modified first-order theory is no 
MOUS T better, since it grossly overestimates the cone loading and gives 
the same general normal force distribution shape. Hybrid 
theory gives results which are very similar to first-order, as would 
—— | be expected, since it is a combination of second-order axial flow 
nd first-order cross flow. On the other hand, modified hybrid 


ons 


theory yields good results for the cone and fair results for the 
frustum 

Second-order shock-expansion theory was used to calculate 
the normal force distribution for all five Mach Numbers. Fol 
wing the technique suggested by Lavender and Deep,* these 
culations were carried out, with the use of a desk calculator, 
inapproximately one hour for each Mach Number. This theory 
isespecially suited to conditions between the ranges of applica 
lity of the perturbation theories and the generalized shock 
expansion theory. It is apparent from Fig. 1 that second-order 
shock-expansion theory gives excellent results for normal force 
listribution. Although total normal force and moment, Figs 
Zand 3, are slightly underestimated, the center of pressure cal 
ulations, Fig. 4, vield accurate results because the predicted 
ormal-force distribution is similar to that obtained from experi 


mental data 
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| tyra’ PICTURE of spotwise transition from laminar to turbu 
lent flow in a boundary layer has been shown to be correct 
in its fundamentals by Schubauer and Klebanoff.2. The inter 
mittency factor y measured by them is a measure of the prob 
ability f(P) that the motion is turbulent at P, but Emmons’ 
assumptions about the source-rate density function g (x, y, ¢) 
do not explain the measured y distributions. Nor is the signifi 
cance of the Gaussian integral curve, fitted by Schubauer and 
Klebanoff to their data, clear. This note directs attention to 
some considerations which suggest an explanation 


It has been shown!’ * that for steady two-dimensional flow 


Wx) = 


l exp g( Xo )dxod yodto 
JR 


the domain of integration being the dependence volume R. In 
zero pressure gradient flow over a flat plate, the experiments of 
Schubauer and Klebanoff have shown that it is justifiable to as 
sume R to be a true cone with straight generators. With this 
assumption, it can be easily shown that the expression for y be 


comes 


¥(x) = xo)*dxo (1) 


°x 
l exp| —(¢/U) } 2(Xo)(x 


where o is Emmons’ spot propagation parameter. One can, in 


principle, obtain g(x) from measured y(x); for successive differ 


entiation of Eq. (1) gives 


g(x) = (U/20)[d°G(x)/dx?| 


where G(x) = In (1 y). In actual practice, however, the 
third derivative involved cannot be obtained with any accuracy 
from a set of experimental points. But it is mteresting to note 
that measured points for G(x) lie strikingly on a parabola with 
vertex at a point x, coinciding with the beginning of transition 
This indicates that the third derivative (or the source rate) is 
nearly zero everywhere except possibly near x;, where it may have 
a singularity 

Schubauer and Klebanoff have also noted that the region of 
laminar breakdown seems limited in extent. This seems reason 
able since the amplified Tollmien-Schlichting waves (from which 
breakdowns presumably occur in some way) might be expected 
to reach critical conditions at a fairly definite x because of the two 
dimensionality of the flow. If therefore, we assume, as a first 
approximation, that all the point-like breakdowns occur (ran 
domly ) across a single line x = x, or in other words that g(x) can 
be approximated by Dirac’s delta function, we can derive quite 
easily that 


v(x) = 1 exp [—(x — x,)*na/U] 


where n is the number of sources produced per unit length and 


unit time across x = x,. If we use the distance \ between the 
points at which y = 0.25 and y = 0.75 (for instance) to charac 
terize the extent of the transition region and put — = (x X)/d 
where # is the point at which y = 0.5, we can express y as a 


unique function of &: 


y¥=1 exp [—0.412 (¢ + 1.3)? (2 
Eq. (2) has been plotted in Fig. 1 along with experimental points 
The agreement is seen to be quite good. It is also seen that if we 


plot F(y) = V} In (1 y){ against — we should get a straight 


# This research was conducted under the guidance of Professor S. Dhawan 
The author is grateful to the National Institute of Sciences of India for the 
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line. Fig. 2 confirms this. The intermittency distribytigy 
against £ is independent of the spot propagation constants an¢ 
explains at once both the similarity and skewness observed j 
experimental y distributions. 

It should be added that investigation of several forms for a(x 
(including normal distributions about x,) indicated that th 
agreement between experiment and calculation was better as ¢ 
g(x) distribution approached a delta function. It is also inter 
esting to note that the virtual origin of the turbulent bounda: 
layer after transition is very near x,. The fully developed turby 
lent boundary layer is a random mix-up of the spots and behave 
like a two-dimensional boundary layer. During transition aj 
the boundary layer is a similar mix-up of spots but for smalle 
times as the turbulent spots alternate with laminar flow. Th 
overall effect seems to be roughly equivalent to a flow where the 
boundary layer alternates between the two-dimensional lamina 
and the two-dimensional turbulent. The thickness of the trans; 
tion boundary layer would be equal to that of the usual laminar 
layer or of the turbulent layer starting from x,, whichever is the 
larger. This indeed is found to be the case. 

In conclusion, it would seem that the theory of Emmons to 
gether with the hypothesis of localized breakdown adequately 
describes the increase of intermittency in the transition region of 
a boundary layer 
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On the Simulation of Random Excitations . 


(Continued from page 649) 


per cent time that X(t) /¢ is less than x—i.e., P| X(t)/o 
< x}. Next, the output voltage of the summing cir- 
cuit is integrated. A time history of this integrated 
voltage resembles a flight of steps which, for a suffi- 
ciently low recording speed, approximates a straight 
line. The measured slope of this line is a direct meas- 
urement of the cumulative probability P} X(t) /o < x}. 
This procedure was used to determine P}X(t)/o < x} 
from approximately x = —2 to x = 2. Negative 
values of x were treated by reversing polarities of the 
random variable. Probability values above x = 2 are 
difficult to distinguish from one; similarly, values be- 
low x = —2 are difficult to distinguish from zero. 
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